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Abstract 

Equations arising in General Relativity are usually too complicated to be solved analytically 
and one has to rely on numerical methods to solve sets of coupled partial differential equations. 
Among the possible choices, this paper focuses on a class called spectral methods where, 
typically, the various functions are expanded onto sets of orthogonal polynomials or functions. 
A theoretical introduction on spectral expansion is first given and a particular emphasis is put 
on the fast convergence of the spectral approximation. We present then different approaches 
to solve partial differential equations, first limiting ourselves to the one-dimensional case, with 
one or several domains. Generalization to more dimensions is then discussed. In particular, the 
case of time evolutions is carefully studied and the stability of such evolutions investigated. 
One then turns to results obtained by various groups in the field of General Relativity by 
means of spectral methods. First, works which do not involve explicit time-evolutions axe 
discussed, going from rapidly rotating strange stars to the computation of binary black holes 
initial data. Finally, the evolutions of various systems of astrophysical interest are presented, 
from supernovae core collapse to binary black hole mergers. 
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1 Introduction 



Einstein's equations represent a complicated set of nonlinear partial differential equations for which 
some exact j30) or approximate analytical solutions are known. But these solutions arc not 
always suitable for physically or astrophysically interesting systems, which require an accurate de- 
scription of their relativistic gravitational field without any assumption on the symmetry or with 
the presence of matter fields, for instance. Therefore, many efforts have been undertaken to solve 
Einstein's equations with the help of computers in order to model relativistic astrophysical objects. 
Within this field of numerical relativity, several numerical methods have been experimented and 
a large variety of them are currently being used. Among them, spectral methods are now increas- 
ingly popular and the goal of this review is to give an overview (at the moment it is written or 
updated) of the methods themselves, the groups using them and the obtained results. Although 
some theoretical framework of spectral methods is given in Sections [2] to HI more details about 
spectral methods can be found in the books by Gottlieb and Orszag [94] , Canuto et al. [Ml ISll EH] , 
Fornberg '75', Boyd f35^ and Hesthaven et al. [f 17| . While these references have of course been 
used for writing this review, they can also help the interested reader to get deeper understanding of 
the subject. This review is organized as follows: hereafter in the introduction, we briefly introduce 
the spectral methods, their usage in computational physics and give a simple example. Section [2] 
gives important notions concerning polynomial interpolation and the solution of ordinary differ- 
ential equations (ODE) with spectral methods. Multi-domain approach is also introduced there, 
whereas some of the multi-dimensional techniques are described in Section [3j The cases of time- 
dependent partial differential equations (PDE), are treated in Section (4] The last two sections are 
then reviewing results obtained using spectral methods: on stationary configurations and initial 
data (Section [5]), and on the time-evolution (Section [6]) of stars, gravitational waves and black 
holes. 



1.1 About spectral methods 

When doing simulations and solving PDEs, one faces the problem of representing and deriving 
functions on a computer, which deals only with (finite) integers. Let us take a simple example of 
a function / : [— f , f] ^ M. The most straightforward way to approximate its derivative is through 

finite-differences methods: first one must setup a grid 

{a;J..=o...JV C [-1,1] 

of A'^ -t- 1 points in the interval, and represent / by its A^ -f 1 values on these grid points 

{fi — fi-'^i)}i=a...N ■ 

Then, the (approximate) representation of the derivative /' shall be, for instance 

Vz < N, f[ = r{xi) Il±l^. (1) 

If we suppose an equidistant grid, so that Vi < A^, x^+i ^ Xi — Ax — 1/N, the error in the 
approximation ^ will decay as Aa; (first-order scheme). One can imagine higher-order schemes, 
with more points involved for the computation of each derivative and, for a scheme of order n, the 
accuracy can vary as (Aa;)" = 1/N^\ 

Spectral methods represent an alternate way: the function / is no longer represented through 
its values on a finite number of grid points, but using its coefficients (coordinates) {cij^^Q jy in a 
finite basis of known functions {^ilj^g n 

N 

/(x)c.^Q<D,(a:). (2) 

i=0 
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A relatively simple case is, for instance, when f{x) is a periodic function of period 2, and the 
^2i{x) = cos(7ria;), $2i+i sin(7ria;) are trigonometric functions. Equation ^ is then nothing but 
the truncated Fourier decomposition of /. In general, derivatives can be computed from the 
Ci's, with the knowledge of the expression for each derivative $^(a;) as a function of {<i>i}j^Q ^■ 
The decomposition ([2|) is approximate in the sense that {<i>j}j^o represent a complete basis 
of some finite-dimensional functional space, whereas / usually belongs to some other infinite- 
dimensional space. Moreover, the coefRcients Ci are computed with finite accuracy. Among the 
major advantages of using spectral methods is the rapid decay of the error (faster than any power 
of l/N, and in practice often exponential e~^), for well-behaved functions (see Section [2.4.4p : one 
therefore has an infinite-order scheme. 

In a more formal and mathematical way, it is useful to work within the methods of weighted 
residuals (MWR, see also Section [2T5|) . Let us consider the PDE 

Lu{x) = s{x) xeU CR"^, (3) 
Buix) = xedU, (4) 

where L is a linear operator, B the operator defining the boundary conditions and s is a source 
term. A function u is said to be a numerical solution of this PDE if it satisfies the boundary 
conditions (|4]) and makes "small" the residual 

R = Lu- s. (5) 

If the solution is searched in a finite-dimensional subspace of some given Hilbert space (any relevant 
space) in terms of the expansion ([2]), then the functions {^i{x)}^^Q are called trial functions 
and, in addition the choice of a set of test functions {^i{x)}^^Q ^ defines the notion of smallness 
for the residual by means of the Hilbert space scalar product 

Vz = 0...iV, (e^,i?) = 0. (6) 

Within this framework, various numerical methods can be classified according to the choice of the 
trial functions: 

• Finite differences: the trial functions are overlapping local polynomials of fixed order 
(lower than iV), 

• Finite elements: the trial fimctions are local smooth functions which are non-zero only on 
sub-domains of C/, 

• Spectral methods: the trial functions are global smooth functions on U . 

Various choices of the test functions define different types of spectral methods, as detailed in 
Section r2.5l Usual choices for the trial functions are (truncated) Fourier series, spherical harmonics 
or orthogonal families of polynomials. 

1.2 Spectral methods in physics 

We do not give here all the fields of physics where spectral methods are being employed, but 
we sketch the variety of equations and physical models that have been simulated with such tech- 
niques. Spectral methods originally appeared in numerical fluid dynamics, where large spectral 
hydrodynamic codes have been regularly used to study turbulence and transition to the turbulence, 
since the seventies. For fully resolved, direct numerical calculations of Navier-Stokes equations, 
spectral methods were often preferred for their high accuracy. Historically, they also allowed for 
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two- or three-dimensional simulations of fluid flows, because of their reasonable computer mem- 
ory requirements. Many applications of spectral methods in fluid dynamics have been discussed 
by Canuto et al. [5S1 [SS], and the techniques developed in that field can be of some interest for 
numerical relativity. 

From pure fluid-dynamics simulations, spectral methods have rapidly been used in connected 
fields of research: geophysics jl90) . meteorology and climate modeling [217j . In this last domain of 
research, they provide global circulation models that are then used as boundary conditions to more 
specific (lower-scale) models, with improved micro-physics. In this way, spectral methods are only 
a part of the global numerical model, combined with other techniques to bring the highest accuracy, 
for a given computational power. Solution to the Maxwell equations can, of course, also be obtained 
with spectral methods and therefore, magneto-hydrodynamics (MHD) have been studied with these 
techniques (see e.g. Hollerbach jll9) ). This has been the case in astrophysics too, where for example 
spectral three-dimensional numerical models of solar magnetic dynamo action realized by turbulent 
convection have been computed [52]. Still in astrophysics, the Kompaneet's equation, describing 
the evolution of photon distribution function in a bath of plasma at thermal equilibrium within the 
Fokker-Planck approximation, has been solved using spectral methods to model the X-ray emission 
of Her X-1 [33j[40]. In the simulations of cosmological structure formation or galaxy evolution, 
many N-body codes rely on a spectral solver for the computation of the gravitational force by the 
so-called particle-mesh algorithm. The mass corresponding to each particle is decomposed onto 
neighboring grid points, thus defining a density field. The Poisson equation giving the Newtonian 
gravitational potential is then usually solved in Fourier space for both fields [118' . 

To our knowledge, the first published results on the numerical solution of Einstein's equations, 
using spectral methods is the spherically-symmetric collapse of a neutron star to a black hole 
by Gourgoulhon in 1991 [W. He used the spectral methods as they have been developed in the 
Meudon group by Bonazzola and Marck l44l . Later studies of fast rotating neutron stars [4lj (sta- 
tionary axisymmetric models), the collapse of a neutron star in tensor-scalar theory of gravity |157j 
(spherically-symmetric dynamical spacetime) and quasi-equilibrium configurations of binary neu- 
tron stars |39 and of black holes [llOJ (three-dimensional and stationary spacetimes) have grown 
in complexity until the three-dimensional time-dependent numerical solution of Einstein's equa- 
tions [37]. On the other hand, the first fully three-dimensional evolution of the whole Einstein 
system has been achieved in 2001 by Kidder et al. |127| . where a single black hole was evolved 
until t ~ 600A/ — ISOOAf , using excision techniques. They used spectral methods as developed in 
the Cornell/Caltech group by Kidder et al. |125j and Pfeiffer et al. |172j . Since then, they have 
focused on the evolution of a binary black hole system, which has recently been simulated until 
the merger and the ring-down by Scheel et al. [186] . Other groups (for instance Ansorg et al. [TU] . 
Bartnik and Norton [21,, Frauendiener [8TJ and Tichy |219) ) have also used spectral methods to 
solve Einstein's equations; Sections [5] and [6] are devoted to a more detailed review of all these 
works. 

1.3 A simple example 

Before entering the details of spectral methods in Sections [21 [3] and [4] let us give here their spirit 
with the simple example of the Poisson equation in a spherical shell: 

A<p = a, (7) 

where A is the Laplace operator (|93)) expressed in spherical coordinates (r, 6, (p) (see also Sec- 
tion [32]) • We want to solve Equation ([7]) in the domain where < i?min 1^ r < i?max, € 
[0,7r], if e [0,27r). This Poisson equation naturally arises in numerical relativity when, for exam- 
ple, solving for initial conditions or the Hamiltonian constraint in the 3-1-1 formalism [97] : the 
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linear part of these equations can be cast into the form ([7]), and the non-hnearities put into the 
source ct, with an iterative scheme on (j). 

First, the angular parts of both fields are decomposed onto a (finite) set of spherical harmonics 
{y/"} (see Section |3X2]): 

l^^^ m=l 

a(r,0,^)c.£ J2 semir)Yr{9,^), (8) 
1=0 m=-e 

with a similar formula relating (jj to the radial functions fem(r). Because spherical harmonics are 
eigenfunctions of the angular part of the Laplace operator, the Poisson equation can be equivalently 
solved as a set of ordinary differential equations for each couple (£, m), in terms of the coordinate 

V(^,to), —71- + — ^ 2 = ShyAr). (9) 

We then map 

r ^ ^ = —^ Zr-^ ' (10) 

-'•'max -'i-miii 

and decompose each field onto a (finite) basis of Chebyshev polynomials ^ (see Sec- 

tion [2X3l) : 

N 



i=0 
N 

flmiO^Y.^'^^MO- (11) 

Each function fim{r) can be regarded as a column- vector Ae„i of its A'^ + 1 coefficients aum in this 
basis; the linear differential operator on the left-hand side of Equation ^ being thus a matrix Lim 
acting on this vector: 

L em Aim — Sim, (12) 

with Stm being the vector of the A'^ -I- 1 coefficients cum of Sim (r) . This matrix can be computed 
from the recurrence relations fulfilled by the Chebyshev polynomials and their derivatives (see 
Section r2. 4. 31 for details). 

The matrix L is singular, because the problem ([7]) is ill-posed. One must indeed specify bound- 
ary conditions at r = i?min and r — -Rmax- For simplicity, let us suppose 

V(0,^), 0(r = i?„in,0,(^) -0(r -ii'^ax,^,'^) = 0. (13) 

To impose these boundary conditions, we shall adopt the tau methods (see Section [2.5.2p : we build 
the matrix L, taking L and replacing the last two lines by the boundary conditions, expressed in 
terms of the coefficients from the properties of Chebyshev polynomials: 

N N 

V(£,m), ^(-l)^a,,„ = ^a,,„ = 0. (14) 

2=0 i=0 

Equations (|14p are equivalent to the boundary conditions p^ . within the considered spectral 
approximation, and they represent the last two lines of L, which can now be inverted and give the 
coefficients of the solution 0. 
If one summarizes the steps: 
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1. Setup an adapted grid for the computatfon of spectral coefRcients (e.g. equidistant in the 
angular directions and Chebyshev-Gauss-Lobatto collocation points, see Section I2.4.3P ; 

2. Get the values of the source a on these grid points; 

3. Perform a spherical-harmonics transform (for example using some available library [152] ). 
followed by the Chebyshev transform (using a Fast Fourier Transform-FFT, or a Gauss- 
Lobatto quadrature) of the source cr; 

4. For each couple of values (I, m), build the corresponding matrix L, with the boundary condi- 
tions and invert the system (using any available linear-algebra package) with the coefficients 
of CT as the right-hand side; 

5. Perform the inverse spectral transform to get the values of </> on the grid points, from its 
coefficients. 

A numerical implementation of this algorithm has been reported by Grandclement et al. [109] , 
who have observed that the error decayed as e~^'"'*^ • e~^, provided that the source a is smooth. 
Machine round-off accuracy can be reached with ^max ^ N ^ 30, which makes the matrix inversions 
of step 4 very cheap in terms of CPU, and the whole method affordable in terms of memory usage. 
These are the main advantages of using spectral methods, as it shall be shown in the following 
sections. 



6 



2 Concepts in One Dimension 



In this section the basic concepts of spectral methods in one spatial dimension are presented. 
Some generalities about the approximation of functions by polynomials are first exposed. The 
basic formulas of the spectral expansion arc then given and two sets of polynomials are discussed 
(Legendre and Chebyshev polynomials). A particular emphasis is put on convergence properties 
(i.e., the way the spectral approximation converges to the real function). 

In Section 12.51 three different methods of solving an ordinary differential equation (ODE) are 
exhibited and applied to a simple problem. Section [2. 61 is concerned with multi-domain techniques. 
After giving some motivations for the use of multi-domain decomposition, four different implemen- 
tations are discussed, as well as their respective merits. One simple example is given, which uses 
only two domains. 

For problems in more than one dimension see Section [3l 

2.1 Best polynomial approximation 

Polynomials are the only functions that a computer can exactly evaluate and so it is natural to 
try to approximate any function by a polynomial. When considering spectral methods, one will 
use global polynomials on a few domains. This is to be contrasted with finite difference schemes, 
for instance, where only local polynomials are considered. 

In this particular section, real functions of [—1, 1] are considered. A theorem due to Weierstrass 
(see for instance [SS]) states that the set P of all polynomials is a dense subspace of all the continuous 
functions on [—1, 1], with the norm IHlg^. This maximum norm is defined as 

ll/lloo= max l/(a;)|. (15) 

This means that, for any continuous function / of [— 1, 1], there exists a sequence of polynomials 
(pi) ,i €N that converges uniformly towards /: 

hm ||/-p,|U,=0. (16) 

This theorem shows that it is probably a good idea to approximate continuous functions by poly- 
nomials. 

Given a continuous function /, the best polynomial approximation of degree iV, is the polyno- 
mial that minimizes the norm of the difference between / and itself: 

11/ -P^lloo = min{||/ -pII^ ,p e Fn} • (17) 

Chebyshev alternate theorem states that for any continuous function /, p*j^ is unique (Theo- 
rem 9.1 of [179] and theorem 23 of 11501 1. There exist iV + 2 points Xi G [-1, 1] such, that the 
error is exactly attained at those points, in an alternate manner: 

/ (x.) - p*^ (x,) = (-1)^+' 11/ - , (18) 

where S = ot S ^ 1. An example of a function and its best polynomial approximation is shown 
on Figure HJ 

2.2 Interpolation on a grid 

A grid X on the interval [—1,1] is a set of iV -|- 1 points Xi E [—1,1], < i < N. These points are 
called the nodes of the grid X . 
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N=2 




Figure 1: Function / — cos"^ {'kx/2) + (x + 1)^ /8 (black curve) and its best approximation of degree 
2 (red curve). The blue arrows denote the 4 points where the maximum error is attained. 



Let us consider a continuous function / and a grid family of grids X with iV + 1 nodes Xi 

X 
N . 

X . 



Then, there exists a unique polynomial of degree iV, /^/, that coincides with / at each node: 



f = f {X^) Q<i<N. (19) 

f is called the interpolant of / through the grid X. f can be expressed in terms of the 
Lagrange cardinal polynomials: 

N 

i=0 

where the £f are the Lagrange cardinal polynomials. By definition, if- is the unique polynomial 
of degree TV, that vanishes at all nodes of the grid X but at Xi, where it is equal to 1. It is easy to 
show that the Lagrange cardinal polynomials can be written as 

N 

^f(-)= n (21) 

Figure [2] shows some examples of Lagrange cardinal polynomials. An example of a function and 
its interpolant on a uniform grid can be seen on Figure [3] 

Thanks to Chebyshev alternate theorem, one can see that the best approximation of degree 
N is an interpolant of the function at A'' + 1 nodes. However, in general, the associated grid 
is not known. The difference between the error made by interpolating on a given grid X can 
be compared to the smallest possible error for the best approximation. One can show that (see 
Prop. 7.1 of [T75j ): 

11/ - I§f\L < (1 + Aiv (X)) 11/ - p^ll^ , (22) 
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Uniform grid N=8 




X 



Figure 2: Lagrange cardinal polynomials (red curve) and l-j on an uniform grid with = 8. 
The black circles denote the nodes of the grid. 



N=4 




Figure 3: Function / = cos^ {'kx/2) + {x + if' /8 (black curve) and its interpolant (red curve)on a 
uniform grid of 5 nodes. The blue circles show the position of the nodes. 
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Figure 4: Function / = ^ (black curve) and its interpolant (red curve) on a uniform grid 

1 + 16x^ 

of 5 nodes (left panel) and 14 nodes (right panel). The blue circles show the position of the nodes, 
where A is the Lebesgue constant of the grid X and is defined as: 

N 

An {X) = max,e[_i,i] ^ (a;) | . (23) 

i=0 

A theorem by Erdos [72] states that, for any choice of grid X, there exists a constant C > 
such that: 

AAr(X) > -ln(A^+l)-C. (24) 

TT 

It immediately follows that A^r oo when N — + oo. This is related to a result from 1914 by 
Faber [73] that states that for any grid, there always exists at least one continuous function / 
which interpolant does not converge uniformly to /. An example of such failure of the convergence 

is show on Figure 31 where the convergence of the interpolant to the function / = ^ is 

1 + IGx"^ 

clearly not uniform (see the behavior near the boundaries of the interval). This is known as the 
Runge phenomenon. 

Moreover, a theorem by Cauchy (Theorem 7.2 of [179J states that, for all functions / e C^^^^\ 
the interpolation error, on a grid AT of iV + 1 nodes is given by 

f (x) - I^f {x) = L_}_).,,x^^ (^) , (25) 

where e E [—1,1]. w^^-^ is the nodal polynomial of X, being the only polynomial of degree + 1, 
with a leading coefhcient 1 and that vanishes on the nodes oi X. It is then easy to show that 

N 

^N+i = n - ■ (26) 

1=0 

On Equation (j25p . one has a priori no control on the term involving f^^^ . For a given function, 
it can be rather large and this is indeed the case for the function / shown on Figure S] (one can 
check, for instance that, |/^^^ (1)| becomes larger and larger). However, one can hope to minimize 
the interpolation error by choosing a grid such that the nodal polynomial is as small as possible. 



10 




Figure 5: Same as Figure 2] but using a grid based on the zeros of Chebyshev polynomials. The 
Runge phenomenon is no longer present. 



A theorem by Chebyshev states that this choice is unique and is given by a grid which nodes 
are the zeros of the Chebyshev polynomial Tjy^i (see Section 12.31 for more details on Chebyshev 
polynomials). With such a grid, one can achieve 

lh^+ilL = ^' (27) 

which is the smallest possible value (see Equation (18), Section 4.2, Chapter 5 of [122]). So, a 
grid based on nodes of Chebyshev polynomials can be expected to perform better that a stan- 
dard uniform one. This is what can be seen on Figured! which shows the same function, and its 
interpolants, as on Figure [4] but with a Chebyshev grid. Clearly, the Runge phenomenon is no 
longer present. It can be checked, that, for this choice of function /, the uniform convergence of 
the interpolant to the function is recovered. The reason is that decreases faster than 

/^+^/ {N + 1)1 increases. Of course Faber's result implies that this can not be true for all the 
functions. There still must exist some functions for which the interpolant does not converge uni- 
formly to the function itself (it is actually the class of functions that are not absolutely continuous, 
like the Cantor function). 



2.3 Polynomial interpolation 
2.3.1 Orthogonal polynomials 

Spectral methods are often based on the notion of orthogonal polynomials. In order to define 
orthogonality, one has to define the scalar product of two functions on an interval [—1,1]. Let us 
consider a positive function w of [—1, 1] called the measure. The scalar product of / and g with 
respect to this measure is defined as: 

if, 9)^= I f{x)g{x)w{x)dx. (28) 

A basis of Pn is then a set of iV -I- 1 polynomials {pn}„^o n- degree n and the polynomials 

are orthogonal: {pi,Pj)^ = ioi i ^ j. 
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N=4 



N=8 




The projection PnI of a function / on this basis is then 

N 



PNf = E (29) 



ri=0 

where the coefficients of the projection are given by 

/„ = jl2JPllL (30) 

The difference between / and its projection goes to zero when N increases: 

ll/-^w/IL ^0 when iV^oo. (31) 

Figure [H] shows the function / = cos'^ {t:x/2) + (x + 1)"^ /8 and its projection on Chebyshev poly- 
nomials (see Section [2.4.31) . for = 4 and = 8, illustrating the rapid convergence of Pat/ 
to/. 

At first sight, the projection seems to be an interesting means of numerically representing a 
function. However, in practice this is not the case. Indeed, to determine the projection of a 
function, one needs to compute the integrals pop . which requires the evaluation of / at a great 
number of points, making the whole numerical scheme impractical. 

2.3.2 Gaussian quadratures 

The main theorem of Gaussian quadratures (see for instance [57,) states that, given a measure w, 
there exist + 1 positive real numbers Wn and + 1 real numbers a;„ G [—1, 1] such that: 

N 

V/eP2w+5, / f{x)w(x)dx = y^f{x^)wn. (32) 

The Wn are called the weights and the Xn are the collocation points. The integer 5 can take several 
values depending on the exact quadrature considered: 
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• Gauss quadrature: 6=1. 

• Gauss-Radau: (5 and xq = —1. 

• Gauss-Lobatto: S = —1 and xq = xn = 1. 

Gauss quadrature is the best choice because it apphes to polynomials of higher degree but 
Gauss-Lobatto quadrature is often more useful for numerical purposes because the outermost 
collocation points coincide with the boundaries of the interval, making it easier to impose matching 
or boundary conditions. More detailed results and demonstrations about those quadratures can 
be found for instance in [57] . 

2.3.3 Spectral interpolation 

As already stated in 12.3.11 the main drawback of projecting a function onto orthogonal polyno- 
mials comes from the difficulty to compute the integrals (|30p . The idea of spectral methods is 
to approximate the coefficients of the projection by making use of the Gaussian quadratures. By 
doing so, one can define the interpolant of a function / by 

N 

In f ^ ^ fnPn (x) , (33) 

where 

fn = — f (xi) Pn (xi) Wi and = '^pl,{Xi)Wi. (34) 

i=0 1=0 

The /„ exactly coincides with the coefficients /„, if the Gaussian quadrature is applicable for 
computing I30p . that is for all / G PAr+a. So, in general, I^f ^ -P/v/ and the difference between 
the two is called the aliasing error. The advantage of using the /„ is that they are computed by 
estimating / at the + 1 collocation points only. 

One can show that InJ and / coincide at the collocation points: iNf(xi) = f {xi) so that 
In interpolates / on the grid which nodes are the collocation points. Figure [7] shows the function 
/ = cos^ (7r/2) + {x + 1)'^ /8 and its spectral interpolation using Chebyshev polynomials, for = 4 
and N = 6. 

2.3.4 Two equivalent descriptions 

The description of a function / in terms of its spectral interpolation can be given in two different, 
but equivalent spaces: 

• in the configuration space if the function is described by its value at the A^ + 1 collocation 
points f (xi). 

• in the coefficient space if one works with the A^ + 1 coefficients fi. 

There is a bijection between both spaces and the following relations enable us to go from one 
to the other: 

• the coefficients can be computed from the values of / (xi) using Equation (j34p . 

• the values at the collocation points are expressed in terms of the coefficients by making use 
of the definition of the interpolant (|33p : 

N 

f (Xi) = ^ fnPn (xi) . (35) 
n=0 
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N=4 



N=6 




Figure 7: Function / = cos^ (7rx/2) + (x + 1)"^ /8 (black curve) and its interpolant I^f on Cheby- 
shev polynomials (red curve), for = 4 (left panel) and N = 6 (right panel). The collocation 
points are denoted by the blue circles and correspond to Gauss-Lobatto quadrature. 



Depending on the operation one has to perform on a given function, it may be more clever to 
work in one space or the other. For instance, the square root of a function is very easily given 
in the collocation space by \/ f (xi), whereas the derivative can be computed in the coefficient 
space if, and this is generally the case, the derivatives of the basis polynomials are known, by 

N 
n=0 

2.4 Usual polynomials 

2.4.1 Sturm Liouville problems and convergence 

The Sturm-Liouville problems are eigenvalue problems of the form: 

— (pu')' + qu = Xwu, (36) 
on the interval [—1,1]. p q and w are real valued functions such that: 

• p{x) is continuously diffcrcntiable, strictly positive and continuous at a; = ±1. 

• q {x) is continuous, non negative and bounded. 

• w (x) is continuous, non negative and integrable. 

The solutions are then the eigenvalues \i and the eigenfunctions Ui {x) . The eigenfunctions are 
orthogonal with respect to the measure w: 

1 

Um {x) Un {x) w (x) dx = for m ^ n. (37) 

1 

Singular Sturm-Liouville problems are particularly important for spectral methods. A Sturm- 
Liouville problem is singular if and only if the function p vanishes at the boundaries x = ±1. One 
can show that, if the functions of the spectral basis are chosen to be the solutions of a singular 
Sturm-Liouville problem, then the convergence of the function to its interpolant is spectral, that 
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is faster than any power-law of N, N being the order of the expansion (see Section 5.2 of |57]). Let 
us precise that this does not necessarily imply that the convergence is exponential. Convergence 
properties are discussed in more details for Legendre and Chebyshev polynomials in Section [2.4.41 

Conversely, it can be shown that spectral convergence is not ensured when considering solutions 
of regular Sturm-Liouville problems [57] . 

In what follows two usual types of solutions of singular Sturm-Liouville problems are considered: 
Legendre and Chebyshev polynomials. 



2.4.2 Legendre polynomials 

Legendre polynomials P„ are eigenfunctions of the following singular Sturm-Liouville problem: 

((l-a:2)p;j' + n(n+l)P„ = 0. (38) 

In the notations of Equation ()36p . p = 1 ~ , q — 0, w = \ and A„ = —n [n + 1). 

It follows that Legendre polynomials are orthogonal on [—1,1] with respect to the measure 
w {x) = 1. Moreover, the scalar product of two polynomials is given by: 

/•I 2 
J -I 2n+l 

Starting from Pt^ — 1 and Pi = x, the successive polynomials can be computed by the following 
recurrence expression: 

(n+l) P„+i (x) = {2n + l) xPn (x) ~ nP^-i {x) . (40) 

Among the various properties of Legendre polynomials, one can note that i) P„ has the same 
parity as n. ii) P„ is of degree n. iii) P„ (±1) = (±1)". iv) P„ has exactly n zeros on [—1, 1]. The 
first polynomials are shown on Figure [51 

The weights and locations of the collocation points associated with Legendre polynomials de- 
pend on the choice of quadrature. 

2 

• Legendre-Gauss: the nodes of Pn+i and Wi = 5-. 

(i-^f)[^A.+i (^«)] 

• Legendre-Gauss-Radau: xq — —1 and the xi are the nodes of Pj^ + Pw+i. The weights are 
given by — -j^^ ^ and Wi = ^. 



(N + iy [N + iy 

ss-Loba 
2 1 



• Legendre-Gauss-Lobatto: a;o = — l,a;jv = l and Xi are the nodes of P^. The weights are 



Wi 



N{N + 1) [P^(x,)f 



These values have no analytic expression, but they can be computed numerically in an efficient 
way. 

Some elementary operations can easily be performed on the coefficient space. Let us assume 



N 



that a function / is given by its coefficients a„ so that / ~ anPn- Then, the coefficients &„ of 

n=0 

N 

H f = bnPn can be found as a function of the a„, for various operators H. For instance. 



15 
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-1 -0.5 0.5 1 

X 



Figure 8: First Legendre polynomials, from Pq to P5 



• if if is the multiplication by x then: 



• if iJ is the derivative: 

N 



6„ = (2n + l) J2 «P- (42) 

p—n+l,p+n odd 

• if _ff is the second derivative: 

N 

6„ = (n + l/2) [p{p+l)-n{n + l)]ap. (43) 

p— n+2,p+n even 

Those kind of relations enable to represent the action of ii as a matrix acting on the vector of 
the a„, the product being the coefficients of Hf, i.e. the 6„. 

2.4.3 Chebyshev polynomials 

Chebyshev polynomials T„ are eigenfunctions of the following singular Sturm-Liouville problem: 

(yr^r^)' + -=^T„.0. (44) 



In the notations of Equation p6p . p = \/l — x"^, q = 0, w = l/Vl — x"^ and A„ — ~n. 

It follows that Chebyshev polynomials are orthogonal on [—1,1] with respect to the measure 
w — 1/ y/l — x^ and the scalar product of two polynomials is 

{Tn, T™) = / -^^dx = J (1 + <5o,0 S^n. (45) 
j-i V 1 — a;^ ^ 
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Figure 9: First Chebyshev polynomials, from Tq to T5. 



Given that Tq = I and Ti = x, the higher order polynomials can be obtained by making use of the 
recurrence 

Tn+l (x) = 2xTn (x) - Tn-1 (x) . (46) 

This implies the following simple properties, i) T„ has the same parity as n. ii) T„ is of degree n. 
iii) T„ (±1) = (±1)". iv) T„ has exactly n zeros on [—1, 1]. The first polynomials are shown on 
Figure M 

Contrary to the Legendre case, both the weights and positions of the collocation points are 
given by analytic formulas: 

(2i + l)7r TT 
• Chebyshev-Gauss: Xi — cos and Wi 



2N + 2 N + 1 

StTZ TT 

• Chebyshev-Gauss-Radau: Xi = cos . The weights are wn = and Wi 

27r 



2iV + 1 



7rz 7r TT 

• Chebyshev-Gauss-Lobatto: Xi = cos—. The weights are wq ~ ~ and Wi — —. 

As for the Legendre case, the action of various linear operators H can be expressed in the 
coefficient space. This means that the coefficients 6„ of Hf are given as functions of the coefficients 
On of /. For instance, 

• if _ff is the multiplication by x then: 

&n = ^ [(1 + <5o«-i)a„_i + a„+i] (n>l). (47) 
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if H is the derivative: 



bn = 



(l+<5on) 



if H is the second derivative: 



N 

E 

p— n+l,p+7i odd 



pap. 



(48) 



N 



J2 pip 



(l + (5on) ^ 

p— n+2,p+7i even 



2 2\ 



(49) 



2.4.4 Convergence properties 

One of the main advantage of spectral method is the very fast convergence of the interpolant InJ 
to the function /, at least for smooth enough functions. Let us consider a C" function /, one can 
place the following upper bounds on the difference between / and its interpolant In/- 



For Legendre: 



For Chebyshev: 



\iNf - /iil2 < 2i/2 J2 Ik''' 



]\frn-l/2 



Co 



k=0 



L2 



!^jv/-/IIl2 < ^Elk 



(fe) 



(fe) 



fe=0 



(50) 

(51) 
(52) 



The Ci are positive constants. An interesting limit of the above estimates concerns a C°° func- 
tion. One can then see that the difference between / and In/ decays faster than any power of N. 
This is the so-called spectral convergence. Let us precise that this does not necessarily imply that 
the error decays exponentially (think about exp for instance). Exponential convergence 

is achieved only for analytic functions, i.e. functions that are locally given by a convergent power 
series. 

An example of this very fast convergence is shown on Figure 1101 The error clearly decays as an 
exponential, the function being analytic, until the level of 10^^^ of the precision of the computation 
is reached (one is working in double precision in this particular case). Figure [TOl illustrates the fact 
that, with spectral methods, very good accuracy can be reached with only a moderate number of 
coefficients. 

If the function is less regular (i.e. not C°°), the error only decays as a power-law, thus making 
the use of spectral method less appealing. It can be easily seen on the worst possible case: the one 
of a discontinuous function. In this case, the estimates (|50ll52p do not even ensure convergence at 
all. On Figure [Til one shows a step function and its interpolant, for various values of N. One can 
see that the maximum difference between the function and its interpolant is not even going to zero 
when N is increasing. This is known as the Gibbs phenomenon. 

Finally, equation (|52p shows that if to > 0, the interpolant converges uniformly to the function. 
The continuous functions that do not converge uniformly to their interpolant, which existence are 
proved by Faber ^3] (see Section [2T2|) . must belong to the functions. Indeed, for the case to = 0, 
Equation (|52p does not prove convergence (neither do Equations (|50p and ([ST]) '). 
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Figure 10: Maximum difference between / = cos^ {'kx/2) + (.t + 1)^ /8 and its interpolant InJ, as 
a function of N . 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



X 



Figure 11: Step function (black curve) and its interpolant, for various values of N . 
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2.4.5 Trigonometric functions 



A detailed presentation of the theory of Fourier transform is beyond the scope of this work. How- 
ever, there is a close link between the so-called discrete Fourier transform and the spectral inter- 
polation and this is briefly outlined here. More details can be found, for instance, in [57] . 
The Fourier transform Pf of a function / of [0, 27r] is given by: 



oo 



Pf (x) = oq + On cos (nx) + bn sin (nx) . (53) 

n=l n=l 

The Fourier transform is known to converge rather rapidly to the function itself, if / is periodic. 

However, the coefficients an and bn are given by integrals of the form / / (x) cos (nx) dx, that 

Jo 

cannot easily be computed (as it was the case for the projection of a function on orthogonal 
polynomials in Section [2.3. 

The solution to this problem is also very similar to the use of the Gaussian quadratures. Let 
us introduce N + 1 collocation points Xi = 2Tri/{N +1). Then, the discrete Fourier coefficients 
with respect to those points are: 

1 ^ 

fc=i 

2 ^ 

an = — ^f (xk) COS (nxk) (55) 
fc=i 

2 ^ 

bn = —^f{xk)sm{nxk) (56) 
and the interpolant Ij^f is then given by: 

N N 

In f {x) — ai} + a„ cos (nx) + 6„ sin [nx) . (57) 

71—1 n—1 

The approximation made by using the discrete coefhcients in place of the real ones is of the 
same nature as the one made when computing the coefficients of the projection ([30]) by means 
of the Gaussian quadratures. Let us mention that, in the case of a discrete Fourier transform, 
the first and last collocation points lie on the boundaries of the interval, as for a Gauss-Lobatto 
quadrature. As for the polynomial interpolation, the convergence of /at/ to / is spectral for all 
periodic and C°° functions. 



2.4.6 Choice of basis 

For periodic functions of [0, 27r[, the discrete Fourier transform is the natural choice of basis. 
If the considered function has also some symmetries, one can use a subset of the trigonometric 
polynomials. For instance, if the function is i) periodic on [0, 27r[ and is also odd with respect to 
x = TT, then it can be expanded on sines only. If the function is not periodic, then it is natural 
to expand it either on Chebyshev or Legendre polynomials. Using Legendre polynomials can be 
motivated by the fact that the associated measure is very simple w (x) — 1. The multi-domain 
technique presented in Section [2.6.5l is one particular example where such property is required. In 
practice, Legendre and Chebyshev polynomials usually give very similar results. 
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2.5 Spectral methods for ODEs 



2.5.1 Weighted residual method 

Let us consider a differential equation of the following form 

Lu{x) = S{x), a;e[-l,l], (58) 

where L is a linear second-order differential operator. The problem admits a unique solution once 
appropriate boundary conditions are prescribed at a; = 1 and x = —1. Typically, one can specify 
i) the value of u (Dirichlet-type) ii) the value of its derivative dxU (Neumann-type) iii) a linear 
combination of both (Robin- type). 

As for the elementary operations presented in Section [2. 4. 21 and 12. 4. 31 the action of L on w can 
be expressed by a matrix Lij. If the coefficients of u with respect to a given basis are the iii, then 
the coefficients of Lu are 

N 

^L.jUj. (59) 

Usually, the Lij can easily be computed by combining the action of elementary operations like the 
second derivative, the first derivative, the multiplication or division by x (see Sec. I2.4.2l and l2.4.3l 
for some examples). 

A function u is an admissible solution of the problem if and only if i) it fulfills the boundary 
conditions exactly (up to machine accuracy) ii) it makes the residual R = Lu — S small. In the 
weighted residual method, one considers a set of + 1 test functions {6i}„=o n i^^^ The 
smallness of R is enforced by demanding that 

iR,^k)=Oyk<N. (60) 

As N increases, the obtained solution is closer and closer to the real one. Depending on the choice 
of the test functions and the way the boundary conditions are enforced, one gets various solvers. 
Three classical examples are presented below. 

2.5.2 The Tau- method 

In this particular method, the test functions coincide with the basis used for the spectral expansion, 
for instance the Chebyshev polynomials. Let us denote Ui and Si the coefficients of the solution u 
and of the source S. 

Given the expression of Lu in the coefficient space (j59p and the fact that the basis polynomials 
are orthogonal, the residual equations ([60]) are expressed as 

N 

LniUi = Sn, Vn < N, (61) 

i=0 

the unknowns being the Ui. However, as such, this system does not admit a unique solution, due 
to the homogeneous solutions of L (i.e. the matrix associated with L is not invertible) and one 
has to impose the boundary conditions. In the Tau-method, this is done by relaxing the last two 
equations (j6ip (i.e. for n ^ N — 1 and n = N) and by replacing them by the boundary conditions 
at X = —1 and x = I. 

The Tau-method thus ensures that Lu and S have the same coefficients, excepted the last ones. 
If the functions are smooth, then their coefficients should decrease in a spectral manner and so 
the "forgotten" conditions are less and less stringent as N increases, ensuring that the computed 
solution converges rapidly to the real one. 
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Figure 12: Exact solution (|64p of Equation (|62p (blue curve) and the numerical solution (red 
curves) computed by means of a Tau- method, for = 4 (left panel) and N — 8 (right panel). 



As an illustration, let us consider the following equation 

d^u du , , , 4e 

_-4-+4. = exp(.)-^^ (62) 

with the following boundary conditions 

u(a; = -l) = and u(a; = l) = 0. (63) 
The exact solution is analytic and is given by 

/ \ / \ sinh(l) /„ X e 

Figure [T^ shows the exact solution and the numerical one, for two different values of N. One 
can note that the numerical solution converges rapidly to the exact one, the two being almost 
indistinguishable for N as small as iV = 8. The numerical solution exactly fulfills the boundary 
conditions, no matter the value of N. 



2.5.3 The collocation method 

The collocation method is very similar to the Tau-method. They only differ in the choice of test 
functions. Indeed, in the collocation method one uses continuous function that are zero at each 
but one collocation point. They are indeed the Lagrange cardinal polynomials already seen in 
Section 12.21 and can be written as (xj) — dtj. With such test functions, the residual equations 
(pH)) are 

Lu{xn) ^ Sixr,), Vn<iV. (65) 

The value of Lu at each collocation point is easily expressed in terms of u by making use of 
(|59|) and one gets: 

N N 

Li-jUjTi {xn) = S (xn) , Vn < N. (66) 

i=Q j=0 

Let us note that even if the collocation method imposes that Lu and S coincide at each collo- 
cation point, the unknowns of the system written in the form (|66)) are the coefficients Un and not 
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N=4 



N=8 




Figure 13: Same as Figure [T^ but for the coUocation method. 

the u As for the Tau- method, the system (j66p is not invertible and boundary conditions must 
be enforced by additional equations. In this case, the relaxed conditions arc the two associated 
with the outermost points, i.e. n ^ and n = N, which are replaced by appropriate boundary 
conditions to get an invertible system. 

Figure [T51 shows both the exact and numerical solutions for Equation (|62p. 



2.5.4 Galerkin method 

The basic idea of Galerkin method is to seek the solution u as a sum of polynomials Gi that 
individually verify the boundary conditions. Doing so, u automatically fulfills those conditions and 
they do not have to be imposed by additional equations. Such polynomials constitute a Galerkin 
basis of the problem. For practical reasons, it is better to chose a Galerkin basis that can easily 
be expressed in terms of the original orthogonal polynomials. 

For instance, with the boundary conditions (j63p . one can choose: 

G2fc(x) = T2k+2ix)-To{x) (67) 
G2fe+i(x) = T2k+3ix)-Ti{x) (68) 

More generally, the Galerkin basis relates to the usual ones by means of a transformation matrix 

N 

G,^Y^ MjiTj, yi<N-2. (69) 

3=0 

Let us mention that the matrix M is not square. Indeed, to maintain the same degree of approx- 
imation, one can consider only — 1 Galerkin polynomials, due to the two additional conditions 
they have to fulfill (see for instance Equations (|67ll68p ). One can also note that, in general, the d 
are not orthogonal polynomials. 

The solution u is sought in terms of the coefficients iif on the Galerkin basis: 

N-2 

u{x)^Y.^kG'^ix)- (70) 

k=0 
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Figure 14: Same as Figure [T^ but for the Galerkin method. 



By making use of Equations (|59p and ([55]) one can express Lu in terms of the uf: 

N-2 N N 

Lu (a:) = E E E M,kL^,T, (x) . (71) 
k=Q i=o j=a 

The test functions used in the Galerkin method are the Gi themselves, so that the residual 
system reads: 

(Lu,G„) = (5,G„), Vn<iV-2 (72) 

where the left-hand-side is computed by means of (|7ip and by expressing the Gi in terms of the Ti 
by (|69p. Concerning the right-hand-side, the source itself is not expanded on the Galerkin basis, 
given that it does not fulfill the boundary conditions. Putting all the pieces together, the Galerkin 
system reads: 

N-2 N N N 

E E E ^'^^nM.uUj {T,\Ti)^Y. 1^^) ' Vn < TV - 2. (73) 

k=0 1=0 i=0 i=0 

This is a system of iV — 1 equations for the — 1 unknowns uf and it can be directly solved, 
because it is well-posed. Once the uf are known, one can obtain the solution in terms of the usual 
basis by making, once again, use of the transformation matrix: 

AT /N-2 \ 

"(^) = E E^^»"« p- (74) 

j=0 \n=0 / 

The solution obtained by the application of this method to the Equation (|62p is shown on 
Figure [Til 

2.5.5 The methods are optimal 

A numerical method is said to be optimal if it does not introduce an additional error to the one 
that would be done by interpolating the exact solution of a given equation. 

Let us call Uoxact such exact solution, unknown in general. Its interpolant is /ArWoxact and the 
numerical solution of the equation is Unum- The numerical method is then optimal if and only if 
ll-^AfWexact - ""exactHoo ^nd ||unum " "exactHoo behave in the Same manner when iV ^ oo. 
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N 



Figure 15: Difference between the exact solution of Equation (|62p and its interpolant (black 
curve) and between the exact and numerical solutions for i) the tau method (green curve and circle 
symbols) ii) the collocation method (blue curve and square symbols) iii) the Galerkin method (red 
curve and triangle symbols). 

In general, optimality is difficult to check because both Woxact and its interpolant are unknown. 
However, for the test problem proposed in Section 12.5.21 this can be done. Figure [15] shows the 
maximum relative difference between the exact solution (|64p and its interpolant and the various 
numerical solutions. All the curves behave in the same manner as A'^ increases, indicating that the 
three methods previously presented are optimal (at least for this particular case) . 

2.6 Multi-domain techniques for ODEs 
2.6.1 Motivations and setting 

A seen in Section r2. 4. 41 spectral methods are very efficient when dealing with C°° functions. How- 
ever, they lose some of their appeal when dealing with less regular functions, the convergence to 
the exact functions being substantially slower. Nevertheless, the physicist has sometimes to deal 
with such fields. This is the case for the density jump at the surface of strange stars or the for- 
mation of shocks to mention only two examples. In order to maintain spectral convergence, one 
then needs to introduce several computational domains such that the various discontinuities of the 
functions lie at the interface between the domains. Doing so, in each domain, one only deals with 
C°° functions. 

Multi-domain techniques can also be valuable when dealing with a physical space either too 
complicated or too large to be described by single domain. Related to that, one can also use several 
domains to increase the resolution in some parts of the space where more precision is required. 
This can be easily done by using a different number of basis functions in different domains. One 
then talks about fixed-mesh refinement. 

Parallelism can also be a reason why several domains may be used. Indeed, one could set a 
solver, dealing with each domain on a given processor, and interprocessor communication is then 



25 



only used for matching the solution across the various domains. The algorithm of Section [2.6.41 is 
well adapted to such purpose. 

In the following, four different multi-domain methods are presented to solve an equation of 
the type Lu = S* on [— 1, 1]. L is a second order linear operator and S a given source function. 
Appropriate boundary conditions are given at the boundaries a; = — 1 and x — 1. 

For simplicity the physical space is split into two domains: 

• first domain: a; < described by xi = 2a; + 1, xi G [— 1, 1], 

• second domain: a: > described by a;2 = 2a; — 1, a;2 G [— 1, 1]. 

N 

If a; < 0, a function u is described by its interpolant in terms of a;i: Inu (a;) = u\Ti {xi (a;)). 

1=0 

The same is true for a; > with respect to the variable X2- Such setting is obviously appropriate 
to deal with problems where discontinuities occur at a; = 0, that is a:i = 1 and X2 = —1. 

2.6.2 Multi-domain tau method 

As for the standard tau- method (see Section r2. 5. 2p and in each domain, the test functions are the 
basis polynomials and one writes the associated residual equations. For instance in the domain 
a; < one gets: 

N 

(T„,i?)=0^^L„Mi = 4 Vn<iV, (75) 

1=0 

the being the coefficients of the source and Lij the matrix representation of the operator. As 
for the one-domain case, one relaxes the last two equations, keeping only iV — 1 equations. The 
same is done in the second domain. 

Two supplementary equations are enforced to ensure that the boundary conditions are fulfilled. 
Finally, the operator L being of second order, one needs to ensure that the solution and its first 
derivative are continuous at the interface a; = 0. This translates to a set of two additional equations 
involving both domains. 

So, one considers 

• — 1 residual equations in the first domain, 

• A^ — 1 residual equations in the second domain, 

• 2 boundary conditions, 

• 2 matching conditions, 

for a total of 2A^-|-2 equations. The unknowns are the coefficients of u in both domains (i.e. the 
u\ and the uf), that is 2A^ -I- 2 unknowns. The system is well posed and admits a unique solution. 

2.6.3 Multi-domain collocation method 

As for the standard collocation method (see Section I2.5.3P and in each domain, the test functions 
are the Lagrange cardinal polynomials. For instance in the domain a; < one gets: 

N N 
1=0 j=0 
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Lij being the matrix representation of the operator and xin the n collocation point in the first 
domain. As for the one-domain case, one relaxes the two equations corresponding to the boundaries 
of the domain, keeping only A'^ — 1 equations. The same is done in the second domain. 

Two supplementary equations are enforced to ensure that the boundary conditions are fulfilled. 
Finally, the operator L being second-order, one needs to ensure that the solution and its first 
derivative are continuous at the interface x = Q. This translates as a set of two additional equations 
involving the coefhcients in both domains. 

So, one considers 

• — 1 residual equations in the first domain, 

• A^ — f residual equations in the second domain, 

• 2 boundary conditions, 

• 2 matching conditions, 

for a total of 2A^-|-2 equations. The unknowns are the coefficients of u in both domains (i.e. the 
u\ and the uf), that is 2A'^-|- 2 unknowns. The system is well posed and admits a unique solution. 

2.6.4 Method based on the homogeneous solutions 

The method exposed here proceeds in two steps. First, particular solutions are computed in each 
domain. Then, appropriate linear combination with the homogeneous solutions of the operator L 
are performed to ensure continuity and impose boundary conditions. 

In order to compute particular solutions, one can rely on any of the methods exposed in 
Section 12.51 The boundary conditions at the boundary of each domain can be chosen (almost) 
arbitrarily. For instance one can use, in each domain, a collocation method to solve Lu = 5, 
demanding that the particular solution Upart is zero at both end of each intervals. 

Then, in order to have a solution in the whole space, one needs to add homogeneous solutions 
to the particular ones. In general, the operator L is second-order and it admits two independent 
homogeneous solutions g and /i, in each domain. Let us note that, in some cases, additional 
regularity conditions can reduce the number of available homogeneous solutions. The homogeneous 
solutions can either be computed analytically if the operator L is simple enough or numerically 
but one then needs to have a method for solving Lu — 0. 

In each domain, the physical solution is a combination of the particular solution and the 
homogeneous ones of the type: 

u = Upart + ag + I3h, (77) 

where a and f3 are constants that must be determined. In the two domains case, we are left with 4 
unknowns. The system they must verify is composed of i) 2 equations for the boundary conditions 
ii) 2 equations for the matching of u and its first derivative across the boundary between the 
two domains. The obtained system is called the matching system and generally admits a unique 
solution. 

2.6.5 Variational method 

Contrary to previously presented methods, the variational one is only applicable with Legendre 
polynomials. Indeed, the method requires that the measure be w (x) = 1. It is also useful to 
extract the second-order term of the operator L and to rewrite it like Lu = u" + H, H being of 
first-order only. 
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In each domain, one writes the residual equation expHcitly: 

(C, i?) = =^ J ^u"dx + J £_ {Hu) da; = y ^Sdx. (78) 

The term involving the second derivative of u is then integrated by parts: 

[^u'] ~ j ^'u'dx + y e {Hu) da; = y ^^dx. (79) 

The test functions are the same as the ones used for the collocation method, i.e. functions 
being zero at all but one collocation point, in both domains (c? — 1,2): {xdj) — Sij. By making 
use of the Gauss quadratures, the various parts of Equation ((79)) can be expressed as {d — 1,2 
indicates the domain): 

. N N N 

Ci^'dx = '^^nixdi)u {xdi)Wi ='^'^DiJDinW^u{xdJ) (80) 
i=0 i=0 j=0 

. N N 

£_n{Hu)dx = '^(.n{xdi){Hu){xdi)w.i^Wn'^HmU{xdi) (81) 
i=0 i=0 
N 

/ ^nSdx = ^ Cn (Xdi) S {Xdi) Wi^ S [xdn) Wn, (82) 

where Dij (resp. Hij) represent the action of the derivative (resp. of H) in the configuration space 

N 

g' {xdk) = ^Dkjg{xdj) (83) 

N 

{Hg){xdk) = ^i/fc,g(x,,). (84) 
i=o 

For points strictly inside each domain, the integrated term [^u'] of Equation (|79p vanishes and 
one gets equations of the form: 

N N N 
- ^ ^ DijDinWiU [Xdj) + W„ ^ HniU {xdi) = 5" {xdn) Wn- (85) 
1=0 j"=0 i=0 

This is a set of TV — 1 equations for each domains {d = 1,2). In the above form, the unknowns are 
the u {xdi), i.e. the solution is sought in the configuration space. 

As usual, two additional equations are provided by appropriate boundary conditions at both 
ends of the global domain. One also gets an additional condition by matching the solution across 
the boundary between the two domains. 

The last equation of the system is the matching of the first derivative of the solution. However, 
instead of writing it "explicitly" , this is done by making use of the integrated term in Equation ([79]) 
and this is actually the crucial step of the whole method. Applying Equation ([75]) to the last point 
xiN of the first domain, one gets: 

N N N 

u' (a^i = 1) = ^ ^ D^jDmrnu (xij) -WN^ HniU (xu) + S (xin) wn- (86) 

i=0 j=0 i=0 
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The same can be done with the first point of the second domain, to get u' {x2 = —1) and the last 
equation of the system is obtained by demanding that u' {xi — \) — u' {x2 = —1) and relates the 
values of u in both domains. 

Before finishing with the variational method, it may be worthwhile to explain why Legendre 
polynomials are used. Suppose one wants to work with Chebyshev polynomials instead. The 

measure is then w (x) = - When one integrates the term containing u" by part one then 

Vl — x'^ 

gets 

j -u"fwdx = [-u'fw] + j u'f'w'dx (87) 

Because the measure is divergent at the boundaries, it is difficult, if not impossible, to isolate the 
term in u' . On the other hand, this is precisely the term that is needed to impose the appropriate 
matching of the solution. 

2.6.6 Merits of the various methods 

From a numerical point of view, the method based on an explicit matching using the homogeneous 
solutions is somewhat different from the two others. Indeed, one has to solve several systems in a 
row but each one is of the same size than the number of points in one domain. This splitting of 
the different domains can also be useful for designing parallel codes. On the contrary, for both the 
variational and the tau method one has to solve only one system, but its size is the same as the 
number of points in whole space, which can be quite large for many domains settings. However, 
those two methods do not require to compute the homogeneous solutions, computation that could 
be tricky depending on the operators involved and on the number of dimensions. 

The variational method may seem more difficult to implement and is only applicable with 
Legendre polynomials. However, on the mathematical grounds, it is the only method which is 
demonstrated to be optimal. Moreover, some examples have been found where the others methods 
are not optimal. It remains true that the variational method is very dependent on both the shape 
of the domains and the type of equation that needs to be solved. 

The choice of one method or another thus depends on the particularity of the situation. As 
for the mono-domain space, for simple tests problems, the results are very similar. Figure [16] 
shows the maximum error between the analytical solution and the numerical one for the four 
different methods. All errors decay exponentially and reach machine accuracy with the roughly 
same number of points. 
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Figure 16: Difference between tlie exact and numerical solutions of the following test problem. 
+ 4u = S, with S (x < 0) — I and S (x > 0) = 0. The boundary conditions are u (x = —1) = 

dx^ 

and u{x = 1) = 0. The black curve and circles denote results from the multi-domain Tau method, 
the red curve and squares from the method based on the homogeneous solutions, the blue curve 
and diamonds from the variational one, the green curve and triangles from the collocation method. 
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3 Multi-Dimensional Cases 



In principle, the generalization to more than one dimension is rather straightforward if one uses the 
tensor product. Let us first take an example, with the spectral representation of a scalar function 
f{x,y) defined on the square {x,y) E [—1,1] x 1, 1] in terms of Chebyshev polynomials. One 
simply writes 

M N 

/(x,2/) = ^^a,,T,(x)T,(2/), (88) 

i=0 j=0 

with Ti being the Chebyshev polynomial of degree i. The partial differential operators can also be 
generalized, as being linear operators acting on the space Pm^Pn- Simple, linear Partial Differen- 
tial Equations (PDE) can be solved by one of the methods presented in Section 12.51 (Galerkin, tau 
or collocation), on this MA^-dimensional space. The development ([88]) can of course be generalized 
to any dimension. Some special PDE and spectral basis examples, where the differential equation 
decouples for some of the coordinates, shall be given in Section [321 

3.1 Spatial coordinate systems 

Most of interesting problems in numerical relativity involve non-symmetric systems and require the 
use of a full set of three-dimensional coordinates. We briefly review hereafter several coordinate 
sets (all orthogonal) that have been used in numerical relativity with spectral methods. They are 
described through the line element ds^ of the flat metric in the coordinates we discuss. 

• Cartesian (rectangular) coordinates are of course the simplest and most straightforward 
to implement; the line element reads ds^ — dx^ + dy^ + dz^ . These coordinates are regular 
in all space, with vanishing connection which makes them easy to use, since all differential 
operators have simple expressions and the associated triad is also perfectly regular. They are 
particularly well-adapted to cubic-like domains, see for instance |168| 1172] . and [81] in the 
case of toroidal topology. 

• Circular cylindrical coordinates have a line element ds^ ~ dp^ + d(p'^ + dz'^ and exhibit 
a coordinate singularity on the z-axis (p = 0). The associated triad being also singular for 
p = 0, regular vector or tensor fields have components that are multi- valued (depending on 0) 
on any point of the z-axis. As for the spherical coordinates, this can be handled quite easily 
with spectral methods. This coordinate system can be useful for axisymmetric or rotating 
systems, see [TU] . 

• Spherical (polar) coordinates will be discussed more in details in Section [321 Their line 
element reads ds^ — dr^ dO^+r^ sir? 6 dip^ , showing a coordinate singularity at the origin 
(r = 0) and on the axis for which = 0, tt. They are very interesting in numerical relativity 
for the numerous spherical-like objects under study (stars or black hole horizons) and have 
been mostly implemented for shell-like domains [101 11091 11681 1219) and for spheres including 
the origin HHTDlj. 

• Prolate spheroidal coordinates consist of a system of confocal ellipses and hyperbola, 
describing an (x, 2:)-plane, and an angle p giving the position, as a rotation with respect to the 
focal axis |131| . The line element is ds^ — a? (sinh^ p + sin^ [dp^ + dv^^+a^ sinh^ psir? v dip^ . 
The foci are situated at z = ±a and represent coordinate singularities for /i = and — 0, tt. 
These coordinates have been used in [8] with black hole punctures data at the foci. 

• Bispherical coordinates are obtained by rotation of bipolar coordinates around the focal 
axis, with a line element ds^ = (coshry — cosx)"^ (d?7^ -I- dx^ + siv? xd'P^)- As for prolate 
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spheroidal coordinates, the foci situated at z — ±a (rj ±00, x = 0,7r) and more generally, 
the focal axis exhibits coordinate singularities. Still, the surfaces of constant r] are spheres 
situated in the z > 0(< 0) region for 77 > 0(< 0), respectively. Thus these coordinate are 
very well adapted for the study of binary systems and in particular for excision treatment of 
binary black holes [6] . 

3.1.1 Mappings 

Choosing a smart set of coordinates is not the end of the story. As for finite-elements, one would 
like to be able to cover some complicated geometries, like distorted stars, tori, etc. . . or even to be 
able to cover the whole space. The reason for this last point is that, in numerical relativity, one 
often deals with isolated systems for which boundary conditions are only known at spatial infinity. 
A quite simple choice is to perform a mapping from numerical coordinates to physical coordinates, 
generalizing the change of coordinates to [—1, 1], when using families of orthonormal polynomials 
or to [0, 2tt] for Fourier series. 




Figure 17: Regular deformation of the [—1, 1] x [—1, 1] square. 

An example of how to map the [—1,1] x [—1,1] domain can be taken from Canuto et al. [5S[ . 
and is illustrated on Figure 1171 once known the mappings from the four sides (boundaries) of f2 
to the four sides of fl, one can construct a two-dimensional regular mapping H, which preserves 
orthogonality and simple operators (see Chapter 3.5 of [55]). 

The case where the boundaries of the considered domain are not known at the beginning of the 
computation can also be treated in a spectral way. In the case where this surface corresponds to 
the surface of a neutron star, two approaches have been used. First in Bonazzola et al. [38], the 
star (and therefore the domain) is supposed to be "star-like", meaning that there exists a point 
from which it is possible to reach any point on the surface by straight lines that are all contained 
inside the star. To such a point is associated the origin of a spherical system of coordinates, so 
that it is a spherical domain, which is regularly deformed to coincide with the shape of the star. 
This is done within an iterative scheme, at every step once the position of the surface has been 
determined. Then, another approach has been developed in Ansorg et al. [10], using cylindrical 
coordinates. It is a square in the plane (p, z) which is mapped onto the domain describing the 
interior of the star. This mapping involves an unknown function, which is itself decomposed on a 
basis of Chebyshev polynomials, so that its coefficients are part of the global vector of unknowns 
(as the density and gravitational field coefficients). 

In the case of binary black hole systems, Scheel et al. [189] have developed horizon-tracking 
coordinates using results from control theory. They define a control parameter as the relative drift 
of the black hole position, and they design a feedback control system with the requirement that 
the adjustment they make on the coordinates be sufficiently smooth, so that they do not spoil 
the overall Einstein solver. In addition, they use a dual-coordinate approach, so that they can 
construct a comoving coordinate map which tracks both orbital and radial motion of the black 
holes and allows them to successfully evolve the binary. The evolutions simulated in that article 
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are found to be unstable, when using a single rotating-coordinate frame. As a remark, we note 
here the work by Bonazzola et al. [32], where another option is explored: the so-called stroboscopic 
technique to match between an inner rotating domain and an outer inertial one. 



3.1.2 Spatial compactification 

As stated above, the mappings can also be used to include spatial infinity into the computational 
domain. Such a compactification technique is not tied to spectral methods and has already been 
used with finite-differences methods in numerical relativity by e.g. Pretorius |177j . However, due 
to the relatively lower number of degrees of freedom necessary to describe a spatial domain within 
spectral methods, it is easier within this framework to use some resources to describe spatial 
infinity and its neighborhood. Many choices are possible to do so, either choosing directly a family 
of well-behaved functions on an unbounded interval, for example the Hermite functions (see e.g. 
Section 17.4 in Boyd [48]), or making use of standard polynomial families, but with an adapted 
mapping. A first example within numerical relativity was given by Bonazzola et al. |41j . with the 
simple inverse mapping in spherical coordinates 

xe[-l,l]. (89) 



a{x — 1) ■ 



This inverse mapping for spherical "shells" has also been used by other authors Kidder and 
Finn [125] . Pfeiffer et al. [172[ 1168) . or by Ansorg et al. in cylindrical [10] and spheroidal [8] 
coordinates. Many more elaborated techniques are discussed in Chapter 17 of Boyd [48], but to 
our knowledge, none has been used in numerical relativity yet. Finally, it is important to point 
out that, in general, the simple compactification of spatial infinity is not well- adapted to solv- 
ing hyperbolic PDEs and the above mentioned examples were solving only for elliptic equations 
(initial data, see Section [5]). For instance, the simple wave equation (|126p is not invariant under 
the mapping ([89]) . as it has been shown e.g. by Sommerfeld (see [202) . Section 23. E). Intuitively, 
it is easy to see that when compactifying only spatial coordinates for a wave-like equation, the 
distance between two neighboring grid points becomes larger than the wavelength, which makes 
the wave poorly resolved after a finite time of propagation on the numerical grid. For hyperbolic 
equations, is is therefore usually preferred to impose physically and mathematically well-motivated 
boundary conditions at finite radius (see e.g. Friedrich and Nagy [83] . Rinne [180] or Buchman and 
Sarbach [53]). 



3.1.3 Patching in more than one dimension 

The multi-domain (or multi-patch) technique has been presented in section 12.61 for one spatial 
dimension. In Bonazzola et al. [40] or Grandclement et al. [109] . the three-dimensional spatial 
domains consist of spheres (or star-shaped regions) and spherical shells, across which the solution 
can be matched as in one dimensional problems (only through the radial dependence). In general, 
when performing a matching in two or three spatial dimensions, the reconstruction of the global so- 
lution across all domains might need some more care to clearly write down the matching conditions 
(see e.g. [168] . where overlapping as well as non-overlapping domains are used at the same time). 
For example in two dimensions, one of the problems that might arise is the counting of matching 
conditions for corners of rectangular domains, when such a corner is shared among more than three 
domains. In the case of a PDF where matching conditions must be imposed on the value of the 
solution, as well as on its normal derivative (Poisson or wave equation), it is sufficient to impose 
continuity of either normal derivative at the corner, the jump in the other normal derivative being 
spectrally small (see Chap. 13 of Canuto et al. [5S]). 

A now typical problem in numerical relativity is the study of binary systems (see also Sec- 
tions [53] and [63]) for which two sets of spherical shells have been used by Gourgoulhon et al. [100] . 
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Figure 18: Two sets of spherical domains describing a binary neutron star or black hole system. 
Each set is surrounded by a compactified domain of the type (|89l) , which is not displayed 



as displayed on Figure [TH Different approaches have been proposed by Kidder et al. |128j , and 
used by Pfeiffer |168j and Scheel et al. |189j where spherical shells and rectangular boxes are com- 
bined together to form a grid adapted to binary black hole study. Even more sophisticated setups 
to model fluid flows in complicated tubes can be found in [144] . 

Multiple domains can thus be used to adapt the numerical grid to the interesting part (manifold) 
of the coordinate space; they can be seen as a technique close to the spectral element method jl67) . 
Moreover, it is also a way to increase spatial resolution in some parts of the computational domain 
where one expects strong gradients to occur: adding a small domain with many degrees of freedom 
is the analog of fixed-mesh refinement for finite-differences. 

3.2 Spherical coordinates and harmonics 

Spherical coordinates (see Figure [T9)) are well- adapted for the study of many problems in numerical 
relativity. Those include the numerical modeling of isolated astrophysical single objects, like a 
neutron star or a black hole. Indeed, stars' surfaces have spherical-like shapes and black hole 
horizons have this topology too, which is best described in spherical coordinates (eventually through 
a mapping, see Section I3.1.ip . As these are isolated systems in General Relativity, the good 
boundary conditions are imposed at infinity, requiring a compactification of space, which is here 
achieved with the compactification of the radial coordinate r only. 

When the numerical grid does not extend to infinity, e.g. when solving for a hyperbolic PDF, 
the boundary defined by r = const is a smooth surface, on which boundary conditions are much 
easier to impose. Finally, spherical harmonics, which are strongly linked with these coordinates 
can simplify a lot the solution of Poisson-like or wave-like equations. On the other hand, there 
are some technical problems linked with this set of coordinates, as detailed hereafter, but spectral 
methods can handle them in a very efficient way. 

3.2.1 Coordinate singularities 

The transformation from spherical (r,9,ip) to Cartesian coordinates {x,y,z) is obtained by 

X = r sin9 cos (fi, (90) 
y = rsinOsmip, (91) 
z = rcos6'. (92) 

One immediately sees that the origin r — 4=> x = y = z = Ois singular in spherical coordinates 
because neither 9 nor (p can be uniquely defined. The same happens for the z— axis, where 9 = or 
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Figure 19: Definition of spherical coordinates (r, 9, ip) of a point M and associated triad (eV, 60,6*1^), 
with respect to the Cartesian ones. 

TT, and if cannot be defined. Among the consequences is the singularity of some usual differential 
operators, as for instance the Laplace operator 

Here, the divisions by r at the center, or by sm9 on the z-axis look singular. On the other hand, 
the Laplace operator, expressed in Cartesian coordinates is a perfectly regular one and, if it is 
applied to a regular function, it should give a well-defined result. So the same should be true if one 
uses spherical coordinates: the operator ([M]) applied to a regular function should yield a regular 
result. This means that a regular function of spherical coordinates must have a particular behavior 
at the origin and on the axis, so that the divisions by r or sinO appearing in regular operators are 
always well-defined. If one considers an analytic function of the (regular) Cartesian coordinates 
f{x, y, z), it can be expanded as a series of powers of x, y and z, near the origin 

/(x,2/,z)= ^ a„p,xV2'- (94) 

Replacing the coordinate definitions ([TOjl - iP^ into this expression gives 

/(r, e,ip)=J2 anpgr"+P+« cos" e sin"+f 9 cos" ip sin^ ip; (95) 

n,p,q 

and rearranging the terms in ip: 

f{r,6,ip)= 6™p,rl"l+2p+9sinl"l+2f 6'cos«6'e""'^. (96) 

With some transformations of trigonometric functions in 0, one can express the angular part in 
terms of spherical harmonics Y^{9, (p), see Section [3.2.2[ with £ — \m\ + 2p + q and obtain the two 
following regularity conditions, for a given couple {£, m): 
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• near 9 = 0, a. regular scalar field is equivalent to f{9) ~ sin'™' 0, 



near r = 0, a regular scalar field is equivalent to /(r) 



r 



In addition, the r-dependence translates into a Taylor series near the origin, with the same parity 
as i . More details in the case of polar (2D) coordinates are given in Chapter 18 of Boyd [48] , 

If we go back to the evaluation of the Laplace operator , it is now clear that the result is 
always regular, at least for £ > 2 and m > 2. We detail the cases of ^ = and £ = 1, using the 
fact that spherical harmonics are eigenfunctions of the angular part of the Laplace operator (see 
Equation (|103p V For i — the scalar field / is reduced to a Taylor series of only even powers of 
r, therefore the first derivative contains only odd powers and can be safely divided by r. Once 
decomposed on spherical harmonics, the angular part of the Laplace operator (|93p acting on the 
£ = 1 component reads — 2/r^, which is a problem only for the first term of the Taylor expansion. 
On the other hand, this term cancels with the f providing a regular result. This is the general 
behavior of many differential operators in spherical coordinates; when applied to a regular field, 
the full operator gives a regular result, but single terms of this operator, may give singular results 
when computed separately; these singularities canceling between two different terms. 

As this may seem an argument against the use of spherical coordinates, let us stress that 
spectral methods are very powerful in evaluating such operators, keeping everything finite. As an 
example, we use Chebyshev polynomials in ^ for the expansion of the field /(r = a^), a being a 
positive constant. From the recurrence relation on Chebyshev polynomials (j46p . one has 



V„>0, Z^i^ = 2T:„(e)-^^^, (97) 
which recursively gives the coefficients of 

form those of /(^). The computation of this finite part g{^) is always a regular and linear operation 
on the vector of coefficients. Thus, the singular terms of a regular operator are never computed, 
but the result is the good one, as if the cancellation of such terms had occurred. Moreover, from 
the parity conditions it is possible to use only even or odd Chebyshev polynomials, which simplifies 
the expressions and saves computer time and memory. Of course, relations similar to Equation (|97p 
exist for other families of orthonormal polynomials, as well as relations to divide by sm9 a function 
developed on a Fourier basis. The combination of spectral methods and spherical coordinates is 
thus a powerful tool for accurately describing regular fields and differential operators inside a 
sphere [44]. To our knowledge, this is the first reference showing that it is possible to solve PDEs 
with spectral methods inside a sphere, including the three-dimensional coordinate singularity at 
the origin. 

3.2.2 Spherical harmonics 

Spherical harmonics are the pure angular functions 



Yri9,^) = f-^^^^Pricos9)e^--, (99) 

where £>0 and \m\ < i. P™[cos9) are the associated Legendre functions defined by 

_ (^ + m)! 1 e 
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for m > 0. The relation 

pn^) - ^^pn^) (101) 

gives the associated Legendre functions for negative m; note that the normahzation factors can 
vary in the hterature. This family of functions have two very important properties. First, they 
represent an orthogonal set of regular functions defined on the sphere; thus any regular scalar field 
if) defined on the sphere can be decomposed onto spherical harmonics 

+OC rn—i 

f{o,v)^J2 E fimYne,^). (102) 

i=0 m=-e 

Since they are regular, they automatically take care of the coordinate singularity on the z-axis. 
Then, they are eigenfunctions of the angular part of the Laplace operator (noted here Ag^): 

V(.,H A..>^r(^,.):=^ + ^^ + -^^ = -.(.+ l)yr(^,.), (103) 

the associated eigenvalues being —£{£ + 1). 

The first property makes the description of scalar fields on spheres very easy: spherical har- 
monics are used as decomposition basis within spectral methods, for instance in geophysics or 
meteorology, and by some groups in numerical relativity [21] 11091 1219] , However, they could 
be more broadly used in numerical relativity, for example for Cauchy-characteristic evolution or 
matching [2281 115] , where a single coordinate chart on the sphere might help in matching quan- 
tities. They can also help to describe star-like surfaces being defined by r = h{6,ip), as event or 
apparent horizons [1531 123| [2], The search for apparent horizons is also made easier: since the 
function h verifies a two-dimensional Poisson-like equation, the linear part can be solved directly, 
just by dividing by —£{£ + 1) in the coefficient space. 

The second property makes the Poisson equation 

A<^(r,0,^) =a(r,0,^) (104) 

very easy to solve (see Section ll.3p . If the source a and the unknown are decomposed onto 
spherical harmonics, the equation transforms into a set of ordinary differential equations for the 
coefficients (see also 11091): 

. d'^<t)ini , 2d(l)tm £{£+l)(f)e,n ..„f.s 

V(£, m) , „ H 5 = cr£m. (105) 

ar^ r ar r-^ 

Then, any ODE solver can be used for the radial coordinate: spectral methods of course (see 
Section but other ones have been used too (see e,g, Bartnik et al. [5ni[ll]). The same technique 
can be used to advance in time the wave equation with an implicit scheme and Chebyshev-tau 
method for the radial coordinate [HI I158[ . 

The use of spherical harmonics decomposition can be regarded as a basic spectral method, as 
the Fourier decomposition. There are therefore publicly available "spherical harmonics transforms" 
which consist of a Fourier transform in the (/^-direction and a successive Fourier and Legendre 
transform in the 6'-direction. A rather efficient one is the SpharmonicsKit/S2Kit [152] . but writing 
one's own functions is also possible [M] . 



3.2.3 Tensor components 

All the discussion in Sections |3.2.1I|?IT^ has been restricted to scalar fields. For vector, or more 
generally tensor fields in three spatial dimensions, a vector basis (triad) must be specified to ex- 
press the components. At this point, it is very important to stress out that the choice of the basis 



37 



is independent of the choice of coordinates. Therefore, the most straightforward and simple choice, 
even if one is using spherical coordinates, is the Cartesian triad ^e^, = e,^ = = . With 

this basis, from a numerical point of view, all tensor components can be regarded as scalars and 
therefore, a regular tensor can be defined as a tensor field whose components with respect to this 
Cartesian frame are expandable in powers of x^y and z (as in Bardeen and Piran |19|). Manip- 
ulations and solutions of PDEs for such tensor fields in spherical coordinates are generalization 
of the techniques for scalar fields. In particular, when using the multi-domain approach with do- 
mains having different shapes and coordinates, it is much easier to match Cartesian components 
of tensor fields. Examples of use of Cartesian components of tensor fields in numerical relativity 
include the vector Poisson equation [109] or, more generally, the solution of elliptic systems arising 
in numerical relativity [172] . In the case of the evolution of the unconstrained Einstein system, 
the use of Cartesian tensor components is the general option, as it is done by the Caltech/Cornell 
group pTlfTgg] . 

The use of an orthonormal spherical basis ^e^ — -§^,^8 — \-§g,^ip = rsine ^) Figure [TO]) 

requires some more care, as it is outlined hereafter. The interested reader can also find some details 
in the works by Bonazzola et al. [HJISZ]- Nevertheless, there are systems in General Relativity 
where spherical components of tensors can be useful: 

• When doing excision for the simulation of black holes, the boundary conditions on the excised 
sphere for elliptic equations (initial data) may be better formulated in terms of spherical 
components for the shift or the 3-metric [5^ 11041 1123j . In particular, the component which 
is normal to the excised surface is easily identified with the radial component. 

• Still in 3+1 approach, the extraction of gravitational radiation in the wave zone is made 
easier if the perturbation to the metric is expressed in spherical components, because the 
transverse part is then straightforward to obtain |218l . 

Problems arise because of the singular nature of the basis itself, in addition to the spherical 
coordinate singularities. The consequences are first that each component is a multi- valued function 
at the origin r = or on the z-axis, and then that components of a given tensor are not independent 
one from another, meaning that one cannot in general specify each component independently or set 
it to zero, keeping the tensor field regular. As an example, we consider the gradient — VV of 
the scalar field (f> = x, where x is the usual first Cartesian coordinate field. This gradient expressed 
in Cartesian components is a regular vector field = 1, = 0, = 0. The spherical 

components of V read 

y = sin cos (p, 
— cos 6* cos (/5, 

= -simp, (106) 

which are all three multi-defined at the origin, and the last two on the z-axis. In addition, if is 
set to zero, one sees that the resulting vector field is no longer regular: for example the square of 
its norm is multi-defined, which is not the good property for a scalar field. As for the singularities 
of spherical coordinates, these difhcultics can be properly handled with spectral methods, provided 
that the decomposition basis are carefully chosen. 

The other drawback of spherical components is that usual partial differential operators mix the 
components. This is due to the non-vanishing connection coefficients associated with the spherical 
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flat metric [37]. For example, the vector Laplace operator (VjV-'y) reads 



Q^i r 9r \ 89 sm 6 sm 9 dip J 

d^V^ 2dV^ 1 2 dV ^cosO dV'' 

-ITT- + -^r- + — ^9^^ + + 2 . 

ar"^ r ar \ sm c'(/7 sm 6 dip sm 



^«^^^ + — — + 2— 27^— -— 27;)' (109) 



with Ag^p defined in Equation (|103p . In particular, the r-component (jl07p of the operator involves 
the other two components. This can make the resolution of a vector Poisson equation, which nat- 
urally arises in the initial data problem [5^ of numerical relativity, technically more complicated 
and the technique using scalar spherical harmonics (Section I3.2.2p is no longer valid. One possi- 
bility can be to use vector, and more generally tensor |146[ 12391 12181 IST] spherical harmonics as 
decomposition basis. Another technique might be to build from the spherical components regular 
scalar fields, which can have a similar physical relevance to the problem. In the vector case, one 
can think of the following expressions 

e^W,V\ x^nV\ M = ^'e^J/cV^y^ (110) 

where r — rer denotes the position vector and eijk the third-rank fully antisymmetric tensor. 
These scalars are the divergence, r-component and curl of the vector field. The reader can verify 
that a Poisson equation for transforms into three equations for these scalars, expandable onto 
scalar spherical harmonics. The reason why these fields may be more interesting than Cartesian 
components is that they can have more physical or geometrical meaning. 



3.3 Going further 

The development of spectral methods linked with the problems arising in the field of numerical 
relativity has always been active and still is now. Among the various directions of research one 
can foresee, quite interesting ones might be the beginning of higher-dimensional studies and the 
development of better-adapted mappings and domains, within the spirit of going from pure spectral 
methods to spectral elements \167, 29 . 



3.3.1 More than three spatial dimensions 

There have been some interest for the numerical study of black holes in higher dimensions: as 
well with compactified extra-dimensions [203] , as in brane world models [2001 1132] ; recently, some 
simulations of the head-on collision of two black holes have already been undertaken [230] . With 
the relatively low number of degrees of freedom per dimension needed, spectral methods should 
be very efficient in simulations involving four spatial dimensions, or more. We give here starting 
points to implement 4-dimensional (as needed by e.g. brane world models) spatial representation 
with spectral methods. The simplest approach is to take Cartesian coordinates (a;,y, z,w), but a 
generalization of spherical coordinates (r, 9, p, ^) is also possible and necessitates less computational 
resources. The additional angle ^ is defined in [0, tt], with the following relations with Cartesian 
coordinates 

X = r sin cos sin^, 
y — r sin sin 93 sin ^, 
z — r cos 6* sin ^, 
w — rcos^. 
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The four-dimensional flat Laplace operator appearing in constraint equations [200j reads 



or-^ r or \ o^-' tan ^ 04 sm ^ / 

where Ag^ is the two-dimensional angular Laplace operator (|103p . As in the three-dimensional 
case, it is convenient to use the eigenf unctions of the angular part, which are here 

G^(cosO^'r(cos6')e™'^, (112) 

with TO integers such that |m| < i < k. P™{x) are the associated Legendre functions defined 
by Equation (|100p . G^(x) are the associated Gegenbauer functions 

GiicosO - (sin^e)Gi'^(cosO with Gi'\x) = (113) 

and Gk{x) being the fc-th Gegenbauer polynomial G^^' with A = 1. Since the Gfc are also particular 
case of Jacobi polynomials with a = /3 = 1/2 (see, for example [131) ). Jacobi polynomials are 
also solutions of a singular Sturm-Liouville problem, which ensures fast convergence properties 
(see Section [2.4. ip . The Gk{x) fulfill recurrence relations that make them easy to implement as 
spectral decomposition basis, like the Legendre polynomials. These eigenfunctions are associated 
with the eigenvalues — fc(fc + 2): 

A4 (Gl(cosOPr(cos6i)e""^) = -k{k + 2)Gi{coa^)Pp (cos 6)6'""^ . (114) 

So as in 3+1 dimensions, after decomposing on such a basis, the Poisson equation turns into 
a collection of ODEs in the coordinate r. This type of construction might be generalized to 
even higher dimensions, with the choice of appropriate type of Jacobi polynomials for every new 
introduced angular coordinate. 
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4 Time-Dependent Problems 



From a relativistic point of view, the time coordinate could be treated in the same way as spatial 
coordinates and one should be able to achieve spectral accuracy for the time representation of a 
space-time function f{t, x, y, z) and its derivatives. Unfortunately, this does not seem to be the 
case and, we are neither aware of any efficient algorithm dealing with the time coordinate, nor of 
any published successful code solving any of the PDE coming from the Einstein equations, with the 
recent exception of the 1+1 dimensional study by Hennig and Ansorg [113j . Why is time playing 
such a special role? It is not obvious to find in the literature on spectral methods a complete and 
comprehensive study. A first standard explanation is the difficulty, in general, to predict the exact 
time interval on which one wants to study the time evolution. Then, time discretization errors in 
both finite-differences and spectral methods are typically much smaller than spatial ones. Finally, 
one must keep in mind that, contrary to finite-differences, spectral methods are storing all global 
information about a function on the whole time interval. Therefore, one reason may be that there 
are strong memory and CPU limitations to fully three-dimensional simulations, it is already very 
CPU and memory consuming to describe a complete field depending on 3-1-1 coordinates, even 
with fewer degrees of freedom, as it is the case for spectral methods. But the strongest limitation 
is the fact that, in the full 3-1-1 dimensional case, the matrix representing a differential operator 
would be of very big size; it would therefore be very time-consuming to invert it in a general case, 
even with iterative methods. 

More details on the standard, finite-differences techniques for time discretization are first given 
in Section |4?T] Due to the technical complexity of a general stability analysis, we first restrict the 
discussion of this section to the eigenvalue stability (Section l4.ip . with the following approach: the 
eigenvalues of spatial operator matrices must fall within the stability region of the time-marching 
scheme. Although this condition is only a necessary one and, in general, is not sufficient, it pro- 
vides very useful guidelines for selecting time-integration schemes. A discussion on the imposition 
of boundary conditions in time-dependent problems is given in Section 14.21 Section 14.31 then de- 
tails stability analysis for spatial discretization schemes, with the examples of heat and advection 
equations, before details of a fully-discrete analysis are given for a simple case (Section 14. 4|) . 



4.1 Time discretization 

There have been very few theoretical developments on spectral time discretization, with the ex- 
ception of lerley et al. [121] , were the authors have applied spectral methods in time for the study 
of the Korteweg de Vries and Burger equations, using Fourier series in space and Chebyshev poly- 
nomials for the time coordinate. They observe a time-stepping restriction: they have to employ 
multi-domain and patching techniques (see Section 12. 6p for the time interval, with the size of each 
sub-domain being roughly given by the Courant-Friedrichs-Lewy (CFL) condition. Therefore, the 
most common approach for time representation are finite-differences techniques, which allow for 
the use of many well-established time-marching schemes, and the method of lines (for other meth- 
ods, including fractional stepping, see Fornberg ■79]). Let us write the general form of a first-order 
in time linear PDE: 

Vt>0, Vxe[-l,l], ^^^^ = Lu{x,t), (115) 

where L is a linear operator containing only derivatives with respect to the spatial coordinate x. 
For every value of the time the spectral approximation um{x, t) is a function of only one spatial 
dimension belonging to some finite-dimensional subspace of the suitable Hilbert space 7i, with the 
given spatial norm, associated for example to the scalar product and the weight w introduced 
in Section r2. 3. II Formally, the solution of Equation (IllSp can be written as: 

Va;e[-l,l], u{x,t) = e^^u{x,Q). (116) 



41 



In practice, to integrate time-dependent problems one can use spectral methods to calculate spatial 
derivatives and standard finite-differences schemes to advance in time. 



4.1.1 Method of lines 

At every instant t, one can represent the function u^ix, t) through a finite set UN{t), composed of 
its time-dependent spectral coefficients, or values at the collocation points. We note Ljy the spectral 
approximation to the operator L, together with the boundary conditions, if a tau or collocation 
method is used. Lat is therefore represented as an TV x TV matrix. This is the so-called method of 
lines, which allows one to reduce a PDE to some ODE, after discretization in all dimensions but 
one. The advantage is that many ODE integration schemes are known (Runge-Kutta, symplectic 
integrators, ...) and can be used here. We shall suppose an equally-spaced grid in time, with the 
time-step noted At and = Un{J x At). 

In order to step from to U^~^^, one has essentially two possibilities: explicit and implicit 
schemes. In an explicit scheme, the action of the spatial operator Lff on C^jv lif<jj must be 
computed to explicitly get the new values of the field (either spatial spectral coefficients or values 
at collocation points). A simple example is the forward Euler method: 

U-^+' ^U'^ + AtLNU^, (117) 

which is first-order and for which, as for any explicit schemes, the time-step must is limited by 
the CFL condition. The imposition of boundary conditions is discussed in Section 14.21 With an 
implicit scheme one must solve for a boundary value problem in term of U^~^^ at each time-step: 
it can be performed in the same way as for the solution of the elliptic equation (|62p presented in 
Section [2.5.21 The simplest example is the backward Euler method: 

Un'^' = U^ + AtLNU-^+\ (118) 

which can be re- written as an equation for the unknown U^'^^: 

iI + AtLN)U^+' ^U'^; 

with / the identity operator. Both types of schemes have different stability properties, which can 
be analyzed as follows. Assuming that Ln can be diagonalized in the sense of the definition given 
in (|4.1.3p , the stability study can be reduced to the study of the collection of scalar ODE problems 

^ = X^Uj,, (119) 
where is any of the eigenvalues of Lat in the sense of Equation ()123|) . 



4.1.2 Stability 

The basic definition of stability for an ODE integration scheme is that, if the time-step is lower 
than some threshold, then \\U'^\\ < Ae^'^^*', with the constants A and K independent of the 
time-step. This is perhaps not the most appropriate definition, since in practice one often deals 
with bounded functions and an exponential growth in time would not be acceptable. Therefore, 
an integration scheme is said to be absolutely stable (or asymptotically stable), if HCT^yII remains 
bounded, VJ > 0. This property depends on a particular value of the product x At. For each 
time integration scheme, the region of absolute stability is the set of the complex plane containing 
all the XiAt for which the scheme is absolutely stable. 

Finally, a scheme is said to be A-stable if its region of absolute stability contains the half 
complex plane of numbers with negative real part. It is clear that no explicit scheme can be 



42 




I f" 



ABl 
AB2 
AB3 
AB4 



Real pai't 

Figure 20: Regions of absolute stability for the Adams-Bashforth integration schemes of order 1 
to 4. 



A-stable due to the CFL condition. It has been shown by Dahlquist that there is no linear 
multi-step method of order higher than 2 which is ^-stable. Thus implicit methods are also limited 
in time-step size if more than second-order accurate. In addition, Dahlquist 65] shows that the 
most accurate second-order yl-stable scheme is the trapezoidal one (also called Crank-Nicolson, or 
second-order Adams-Moulton scheme) 



J+i 

N 



At 



NUn) 



(120) 



On Figures [20l and [211 are displayed the absolute stability regions for the Adams-Bashforth and 
Runge-Kutta families of explicit schemes (see for instance [S^). For a given type of spatial linear 
operator, the requirement on the time-step usually comes from the largest (in modulus) eigenvalue 
of the operator. For example, in the case of the advection equation on [—1,1], with a Dirichlet 
boundary condition 



Lu ■ 



du 

Vt, u(l,i) = 0, 



(121) 



and using a Chebyshev-tau method, one has that the largest eigenvalue of Ljv grows in modulus 
as iV^. Therefore, for any of the schemes considered on Figures and the time-step has a 
restriction of the type 

M<0{N-^), (122) 

which can be related to the usual CFL condition by the fact that the minimal distance between 
two points of a (iV-point) Chebyshev grid decreases hke 0(iV~^). Due to the above cited Second 
Dahlquist barrier [66], implicit time marching schemes of order higher than two also have such 
kind of limitation. 
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Figure 21: Regions of absolute stability for the Runge-Kutta integration schemes of order 2 to 5. 
Note that the size of the region is increasing with the order. 



4.1.3 Spectrum of simple spatial operators 

An important issue in determining the absolute stability of a time-marching scheme for the solution 
of a given PDE is the computation of the spectrum (A;) of the discretized spatial operator ()119p . 
As a matter of fact, these eigenvalues are those of the matrix representation of L^r, together with 
the necessary boundary conditions for the problem to be well-posed (e.g. Bnu — 0). If one notes b 
the number of such boundary conditions, each eigenvalue (here, in the case of the tau method) 
is defined by the existence of a non-null set of coefhcients such that 

(Vj) l<j<N-b, {Lnu)j = A,cj, 

Bnu = 0. (123) 

As an example, let us consider the case of the advection equation (first-order spatial derivative) 
with a Dirichlet boundary condition, solved with the Chebyshev-tau method (|121|) . Because of 
the definition of the problem ()123|) . there are TV — 1 "eigenvalues", which can be computed, after 
a small transformation, using any standard linear algebra package. For instance, it is possible, 
making use of the boundary condition, to express the last coefficient as a combination of the other 
ones 

N-l 

CN = -J2 (124) 

One is thus left with a usual eigenvalue problem for a {N — 1) x {N — 1) matrix. Results are 
displayed on Figure [22] for three various values of N. Real parts are all negative: the eigenvalue 
which is not displayed lies on the negative part of the real axis and is much larger in modulus (it 
IS growing like OiN"^)) than the A'' - 1 others. 

This way of determining the spectrum can be, of course, generalized to any linear spatial 
operator, for any spectral basis, as well as to the collocation and Galerkin methods. Intuitively 
from CFL-type limitations, one can see that in the case of the heat equation {Lu — d'^u/dx'^), 
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Figure 22: Eigenvalues of the first derivative-tau operator (|123p for Chebyshev polynomials. The 
largest (in modulus) eigenvalue is not displayed; this one is real, negative and goes as O(iV^). 



explicit time-integration schemes (or any scheme which is not A-stable) shall have a severe time- 
step limitation of the type 

At < OiN-^), (125) 

for both Chebyshev or Legendre decomposition basis. Finally, one can decompose a higher-order 
in time PDE into a first-order system and then use one of the above proposed schemes. In the 
particular case of the wave equation 

it is possible to write a second-order Crank-Nicolson scheme directly jl58j 



Ui-^' = - U-lf' + ^ ( + ) . (127) 



Since this scheme is ^-stable, there is no limitation on the time-step At, but for explicit or higher- 
order schemes this limitation would be At < 0{N^^), as for advection equation. The solution of 
such an implicit scheme is obtained as that of a boundary value problem at each time-step. 

4.1.4 Semi-implicit schemes 

It is sometimes possible to use a combination of implicit and explicit schemes to loosen a time- 
step restriction of the type (|122p . Let us consider as an example the advection equation with 

non-constant velocity on [—1, 1] 

with the relevant boundary conditions, which shall in general depend on the sign of v{x). If on the 
one hand the stability condition for explicit time schemes (|122l) is too strong, and on the other hand 
an implicit scheme is too lengthy to implement or to use (because of the non-constant coefficient 
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v{x)), then it is interesting to consider the semi-implicit two-step method (see also [94) 



..J+l/2 At 7+1/2 _ J At , T-\JJ.J 

'Jn ^^n'^N - '-'n + [-^N - ^n) '^Ni 

U'l^'~^Lt,U-l,^' - + ^ (iA^ - i^^) ^^ (129) 

where and LJj are respectively the spectral approximations to the constant operators —v{l)d/dx 
and —v{—l)d/dx, together with the relevant boundary conditions (if any). This scheme is abso- 
lutely stable if 

^ \ , M - (130) 

^ 7Vmax|u(x)| ^ ' 

With this type of scheme, the propagation of the wave at the boundary of the interval is treated 
implicitly, whereas the scheme is still explicit in the interior. The implementation of the implicit 
part, for which one needs to solve a boundary-value problem, is much easier than for the initial 
operator (|128p because of the presence of only constant-coefficient operators. This technique is 
quite helpful in the case of more severe time-step restrictions (|125[) . for example for a variable 
coefficient heat equation. 



4.2 Imposition of boundary conditions 

The time-dependent PDE (jllSp can be written as a system of ODEs in time either for the time- 
dependent spectral coefficients {ci(t)}j^g jy the unknown function u{x, t) (Galerkin or tau meth- 
ods), or for the time-dependent values at collocation points {u{xi,t)},-^Q ^ (collocation method). 
Implicit time-marching schemes (like the backward Euler scheme (|118p ) are technically very similar 
to a succession of boundary-value problems, as for elliptic equations or Equation ((62|) described 
in Section [^31 The coefficients (or the values at collocation points) are determined at each new 
time-step by the inversion of the matrix of the type / -I- AtL or its higher-order generalization. 
To represent a well-posed problem, this matrix needs in general the incorporation of boundary 
conditions, for tau and collocation methods. Galerkin methods are not so useful if the boundary 
conditions are time-dependent: this would require the construction of a new Galerkin basis at each 
new time-step, which is too complicated and/or time-consuming. We shall therefore discuss in the 
following sections the imposition of boundary conditions for explicit time schemes, with the tau or 
collocation methods. 



4.2.1 Strong enforcement 

The standard technique is to enforce the boundary conditions exactly, i.e. up to machine precision. 
Let us suppose here that the time-dependent PDE (|115p . which we want to solve, is well-posed 
with the boundary condition 

Vi>0, u(a; = l,t) = 5(t), (131) 

where b{t) is a given function. We give here some examples, with the forward Euler scheme (jll7p 
for time discretization. 

In the collocation method, the values of the approximate solution at (Gauss-Lobatto type) 
collocation points {xi}^^^ jy determined by a system of equations: 

yi = O...N-l, U-l+\xi) = U^{x,) + At{LNU^)ix = x,), (132) 
U-j^j+\x = XN = 1} = b{{J+l)At), 

the value at the bomidary {x = 1) is directly set to be the boundary condition. 
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In the tau method, the vector is composed of the + 1 coefBcients {ci( J x At)}^^^ jy 
at the J-th time-step. If we denote by [LnU^) . the i-th coefficient of Lat appHed to U^, then the 
vector of coefficients {cij^^Q ^ is advanced in time through the system: 

yi = O...N -I, Ci{{J+l)xAt) = c,{J X At) + At{LNU^)^ (133) 



CN{{J+l)xAt) = 6((J+l)At)- ^ cfe, 



k=0 

the last equahty ensures the boundary condition in the coefficient space. 
4.2.2 Penalty approach 

As shown in the previous examples, the standard technique consists in neglecting the solution to 
the PDE for one degree of freedom, in configuration or coefficient space, and using this degree of 
freedom in order to impose the boundary condition. However, it is interesting to try and impose 
a linear combination of both the PDE and the boundary condition on this last degree of freedom, 
as it is shown by the next simple example. We consider the simple (time-independent) integration 
over the interval x £ [—1,1]: 

dii 

— = sm{x - 1), and u{l) = 0, (134) 
dx 

where u{x) is the unknown function. Using a standard Chebyshev collocation method (see Sec- 
tion [23T3]), we look for an approximate solution uat as a polynomial of degree N verifying 

Vi = 0...iV-l. '^(xi) = sin(x, -1), 
ax 

^(x.^l) = 0, 
ax 

where the {xi}^^Q are the Chebyshev-Gauss-Lobatto collocation points. 

We now adopt another procedure that takes into account the differential equation at the bound- 
ary as well as the boundary condition, with un verifying (remember that xn — 1): 

V^ = 0...7V-l, ^(x,) = sin(x, -1), (135) 
dx 

'^^(xn) - tun{xn) = sm{xN - 1), 
dx 

where r > is a constant; one notices that taking the limit r +oo, both system become equiv- 
alent. The discrepancy between the numerical and analytical solutions is displayed on Figure [23l 
as a function of that parameter r, when using N = 8. It is clear from that figure that there exists 
a finite value of r (r„iin — 70) for which the error is minimal and, in particular, lower than the 
error obtained by the standard technique. Numerical evidences indicate that Tmin ^ N'^ ■ This is a 
simple example of weakly imposed boundary conditions, with a penalty term added to the system. 
The idea of imposing boundary conditions up to the order of the numerical scheme has been first 
proposed by Funaro and Gottlieb [85] and can be efficiently used for time-dependent problems, as 
illustrated by the following example. For a more detailed description, the interested reader can 
refer to the review article by Hesthaven |115| . 
Let us consider the linear advection equation 

Vxe hl,l],Vt>0, ^ ^ ^ (136) 
Vt>0, u{l,t) = fit), (137) 
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Figure 23: Behavior of the error in the solution of the differential equation (|134p . as a function of 
the parameter r entering the numerical scheme (jl35p . 



where f{t) is a given function. We look for a Legendre collocation method to obtain a solution, and 
define the polynomial Q^{x), which vanishes on the Legendre-Gauss-Lobatto grid points, except 
at the boundary x = 1: 



Q-ix) = 



Thus, the Legendre collocation penalty method uniquely defines a polynomial UN{x,t) through its 
values at Legendre-Gauss-Lobatto collocation points {a;^}^^ 



^O...N 



yi^o...N, 



dt 



du 



N 



dx 



rQ-{x,)iuN{l,t)~ fit)), 



(138) 



where r is a free parameter as in (|135p . For all the grid points, except the boundary one, this 
is the same as the standard Legendre collocation method (V« — . . . N — 1, Q~{xi) — 0). At 
the boundary point x = xn — 1, one has a linear combination of the advection equation and the 
boundary condition. Contrary to the case of the simple integration (|135p . the parameter t here 
cannot be too small: in the limit r ^ 0, the problem is ill posed and the numerical solution diverges. 
On the other hand, we still recover the standard (strong) imposition of boundary conditions when 
r +00. With the requirement that the approximation be asymptotically stable, we get for the 
discrete energy estimate (see details about the technique in Section li . 3 . 21 hereafter ) the requirement 



^ut,{x,,t) — - 

i=0 



- TU%(t, xn)wn < 0. 



Using the property of Gauss-Lobatto quadrature rule (with the Legendre-Gauss-Lobatto weights 
Wi), and after an integration by parts, the stability is obtained if 



r > 



> 



N{N +1) 



(139) 
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It is also possible to treat more complex boundary conditions, as described in Hesthaven and 
Gottlieb |116j in the case of Robin- type boundary conditions (see Section 12.5.11 for a definition) . 
Specific conditions for the penalty coefficient r are derived but the technique is the same: for each 
boundary, a penalty term is added which is proportional to the error on the boundary condition at 
the considered time. Thus, non-linear boundary operators can also be incorporated in an easy way 
(see e.g. the case of the Burgers equation in [115] ). The generalization to multi-domain solution is 
straightforward: each domain is considered as an isolated one, which requires boundary conditions 
at every time-step. The condition is imposed through the penalty term containing the diiference 
between the junction conditions. This approach has very strong links with the variational method 
presented in Section 12.6.51 in the case of time-independent problems. More detailed discussion 
about weak imposition of boundary condition is given in Canuto et al. (Sec. 3.7 of [57] and Sec. 5.3 
of |58j for multi-domain methods). 



4.3 Discretization in space: stability and convergence 

After dealing with temporal discretization, we now turn to another fundamental question of numer- 
ical analysis of initial value problems, which is to find conditions under which the discrete (spatial) 
approximation UN{x,t) converges to the right solution u{x,t) of the PDE (|115[) as iV — > oo and 
t € [0, T]. The time derivative term is treated formally, as could also be treated a source term on 
the right-hand side, that we do not consider here, for better clarity. 

A given spatial scheme to the PDE is said to be convergent if any numerical approximation 
UN{x,t) obtained through this scheme to the solution u{x,t) 

\\Pnu- unWlI N oo. (140) 

Two more concepts are helpful in the convergence analysis of numerical schemes: 

• consistency: an approximation to the PDE II115|) is consistent if Vw G ?i both 



\Pn{Lv - Lnv)\\l2^ 
\Pnv - vn\\l2 — > 



as TV oo; (141) 



• stability: with the formal notations of Equation (I116|) . an approximation to the PDE (IllSp 
is stable if 

I = sup ; 

l^lli?,, 



WN, ||e^"*|| = sup " ,, ,, " " < C{t), (142) 



where C{t) is independent of TV and bounded for t € [0, T]. 



4.3.1 Lax—Richtmyer theorem 

The direct proof of convergence of a given scheme is usually very difficult to obtain. Therefore, a 
natural approach is to use the Lax-Richtmyer equivalence theorem: "a consistent approximation 
to a well-posed linear problem is stable if and only if it is convergent" . Thus, the study of conver- 
gence of discrete approximations can be reduced to the study of their stability, assuming they are 
consistent. Hereafter, we sketch out the proof of this equivalence theorem. 
The time-evolution PDE (jlisp is approximated by 

^ = Lnuj,. (143) 

To show that stability implies convergence, we subtract it to the exact one (|115p 

d{u — upf) 



dt 



— Ln {u ~ un) + Lu — Ljsiu^ 
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and obtain after integration (the dependence on the space coordinate x is skipped) 

u(t)~UAr(t) =e^"*[u(0)-itAr(0)]+ / e^" [Lu(s) - LAru(s)] ds. (144) 

Using the stabihty property (|142p . the norm (L^) of this equation imphes 

||u(i)-'UAr(t)|| <C(t)|lM(0)-itAr(0)|l+ / C{t - s)\\Lu{s) - LNu{s)\\ds. (145) 

JO 

Since the spatial approximation scheme is consistent and C{t) is a bounded function independent 
of N , for a given i e [0, T] the left-hand side goes to zero as — > oo, which proves the convergence. 
Conversely, to show that convergence implies stability, we use the triangle inequality to get 

< |||e^"*u|| - ||e^*w||| < ||e^"*u - e^*w|| 0. 

From the well-posedness lie'^'uH is bounded and therefore ||e^™*u|| is bounded too, independently 
of iV. 

The simplest stability criterion is the von Neumann stability condition: if we define the adjoint 
L* of the operator L, using the inner product, with weight w of the Hilbert space 

y{u,v)eH, {u,Lv)^^ {L*u,v)^, 

then the matrix representation of L* is also the adjoint of the matrix representation of L^. 
The operator Ln is said to be normal if it commutes with its adjoint L^. The von Neumann 
stability condition is that for normal operators, if there exists a constant K independent of N, 
such that 

yi,l<i<N, Re{Xi)<K, (146) 

with (Ai) being the eigenvalues of the matrix Ljy, then the scheme is stable. This condition pro- 
vides an operational technique for checking the stability of normal approximations. Unfortunately, 
spectral approximations using orthogonal polynomials have in general strongly non-normal matri- 
ces Z/jv and therefore, the von Neumann condition cannot be applied. Some exceptions include 
Fourier-based spectral approximations for periodic problems. 

4.3.2 Energy estimates for stability 

The most straightforward technique for establishing the stability of spectral schemes is the energy 
method: it is based on choosing the approximate solution itself as a test function in the evaluation of 
residual (|60p . However, this technique only provides a sufficient condition and in particular, crude 
energy estimates indicating that a spectral scheme might be unstable may be very misleading for 
non- normal evolution operators (see the example in Section 8 of Gottlieb and Orszag _94 ). 

Some sufficient conditions on the spatial operator L and its approximation Ljv are used in the 
literature to obtain energy estimates and stability criteria. Some of them are listed hereafter: 

• If the operator L is semi-bounded: 

37, L + L* < 7/, (147) 

where I is the identity operator; 

• In the parabolic case, if L satisfies the coercivity condition (see also Chap. 6.5 of Canuto et al. 57jl]): 

3A>0,y{u,v), \{Lu,v)\<A\\u\\\\vl (148) 

^Note the difference in sign convention between | 57| and here 
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and the continuity condition: 

3a>0,VM, {Lu,u) < ~a\\u\\'^; (149) 

• In the hyperbolic case, if there exists a constant C > such that 

Vu, ||iu|| < C||ii||, (150) 
and if the operator verifies the negativity condition: 

Vu, {Lu,u)<0. (151) 

As an illustration, we consider hereafter a Galerkin method applied to the solution of Ea. l|115p . 
where the operator L is semi-bounded, following the definition (|147p . The discrete solution is 
such that the residual (|60p estimated on the approximate solution un itself verifies 

'^^"^ Lun,un] =0. (152) 



V dt 

Separating the time-derivative and the spatial operator: 

— hNiml = \iiL + L*)uN{t), UN{t))^ , 
which shows that the "energy" 

lkiv(i)ll' <e''*||ujv(0)||2 (153) 
grows at most exponentially with time. Since uj\i{t) — e^"*ujv(0) for any ujv(0)i we obtain 

||e^"*|| < e5T* (154) 

which gives stability and therefore convergence, provided that the approximation is consistent 
(thanks to Lax-Richtmyer theorem). 

4.3.3 Examples: heat equation and advection equation 

Heat equation 

We first study the linear heat equation 



— - — = 0, with - 1< a; < 1, t > 0, (155) 



du d'^u 
dt dx^ 

homogeneous Dirichlet boundary conditions 

yt > 0, u(-l, t) = u(l, t) = 0, (156) 

and initial condition 

V - 1 < x < 1, u{x, 0) = u°{x). (157) 

In the semi-discrete approach, the Chebyshev collocation method for this problem (see Section [2.5.3|) 
can de devised as follows: the spectral solution UAr(t > 0) is a polynomial of degree N on the inter- 
val [—1, 1], vanishing at the cndpoints. On the other Chebyshev-Gauss-Lobatto collocation points 
{^fc}fe=i N-i (^®^ Section r2.4.3l) . UN{t) is defined through the collocation equations 

Vfc = l...iV-l, ^(xfc,i)-H(xfc,t) = 0, (158) 
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which are thiie ODEs (discussed in Section \4A\ with the initial conditions 

Vfc = 0...7V, UN{xk,0)^u''{xk). (159) 

The stabihty of such a scheme is now discussed, with the computation of an energy bound to the 
sohition. Multiplying the fc-th equation of the system (|158p by ujv(a;fc, t)'Wk^ where {wk}i^^Q pj are 
the discrete weights for the Chebyshev-Gauss-Lobatto quadrature (Section I2.4.3P , and simiming 
over k; one gets: 

1 d ^ ^ d^UN 

■^-^^{uNixk,t)) Wfc - {xk,t)uN{xk,t)wk = 0. (160) 

fe=0 k=0 

Boundary points (fc — 0, N) have been included in the sum since un is zero there from the 
boundary conditions. The product un x d'^UN/dx'^ is a polynomial of degree 27V — 2, so the 
quadrature formula is exact 

^ [Xk,t)UN[Xk,t)Wk = / ^ [Xk,t)UN[Xk,t)w[x)dx; (161) 

k=0 ^ "^^1 ^ 

and integrating by parts twice, one gets the relation 
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w 



Q^2 ixk,t)uNixk,t)w{x)dx ^ J [j'q^ j 'w{x)dx-- J Ujv-^-jdx. (162) 
By the properties of the Chebyshev weight 

and —- = (1 + 2x2) w^, (--^gg^ 



dx^ w \ dx J dx^ 
it is possible to show that 

u%-^-jdx<3j u%w^dx<6j -^^{xk,t)uN{xk,t)w{x)dx\ (164) 



and thus that 

1 ( \ 

-(xfc, i)uAr(xfc, i)w(x)da; > - — I w{x)dx > 0. (165) 

4 J_i \ dx J 



_i dx^ 

Therefore, integrating over the time interval [0,t\ the relation p60p . one obtains the result that 

N N 

'V{uN{xk,t)fwk < V(u°(a;fc))'u'fe <2 max \u'>{x)\\ (166) 

a:e[0.11 

The left-hand side represents the discrete norm of mat (t) 2, but since this is a polynomial of degree 
2N, one cannot apply the Gauss-Lobatto rule. Nevertheless, it has been shown (see e.g. Section 5.3 
of Canuto et al. [57j ) that discrete and L^-norms are uniformly equivalent, therefore: 

t)fw{x) < 2 max \u"{x)f, (167) 



1 



xe[o,i] 



which proves the stability of the Chebyshev-coUocation method for the heat equation. Gonvergence 
can again be deduced from the Lax-Richtmyer theorem, but a detailed analysis (cf. Section 6.5.1 
of Ganuto et al. [57]) shows that the numerical solution obtained by the method described here 
converges to the true solution and one can obtain the convergence rate. If the solution u{x,t) is 
m-times differentiable with respect to the spatial coordinate x (see Section [2.4.4p the energy norm 
of the error decays like iV^~™. In particular, if the solution is C°°, the error decays faster than any 
power of N. 
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Advection equation 

We now study the Legendre-tau approximation to tlie simple advection equation 

Qvb du 

— + — ^Q, with - 1 < x < 1, i > 0, (168) 
at ox 

homogeneous Dirichlet boundary condition 

Vt>0, it(-l,t) = 0, (169) 

and initial condition 

V-l<a;<l, u(x,Q) ^u°{x). (170) 
If we seek the solution as the truncated Legendre series: 

N 

UN{x,t) = y^^ai{t)Pi{x) 

i=0 

by the tau method, then satisfies the equation: 

duN , duN u\u t \ /i7i^ 

Equating coefficients of Pjv on both sides of (|17ip . we get 

dflAT 
"^ = ^- 

Applying the scalar product with u^r to both sides of Equation (|17ip , we obtain 

i|(ll"«f-<.W = -/'.^<— -H(1)<0; 

which implies the following inequality: 

, N-l 



dt 

1=0 



Finally, aN{t) is bounded because it is determined in terms of {ajjj^Q from the boundary 

condition (|169p . and thus stability is proved. In the same way as before, for the heat equation, it 
is possible to derive a bound for the error ||M(a;,t) — ujv(a;,i)||, if the solution u{x,t) is m-times 
differentiable with respect to the spatial coordinate x: the energy norm of the error decays like 
]\ji-m ^ggg g^jgQ Section 6.5.2 of Canuto et al. [57]). In particular, if the solution is C°°, the error 
decays faster than any power of N . 



4.4 Fully-discrete analysis 



Stability issues have been discussed separately for time (Section I4.1.2p and space (Section 14. 3p 
discretizations. The global picture {fully discrete analysis), taking into account both discretizations 
is, in general, very difficult to study. However, it is possible in some particular cases to address 
the problem and in the following lines, we shall perform a fully discrete analysis of the advection 
equation (|168p . when using a Legendre collocation method in space and a forward Euler scheme 
in time. With the notations of Section WA\ the scheme writes 



yx e [-1,1], 



U'^+\x)^U-^{x) 



ox 



dx 



Pn+i{x) 



(173) 
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where the last term imposes the boundary condition VJ, U^{x = —1) = 0. We consider this 
relation at the Legendre-Gauss collocation points {{xi}^^^ which are zeros of Pjv+i(a;); the 
square of this expression taken at these collocation points gives 



V» e [0,iV], {Ui+\xi)f = {ui,{x,)f+Ae 



N 



dx 



dx 



we multiply by (1 — Xi)wi, where the {wi}-^^ jy are the Legendre-Gauss weights, and sum over i 
to obtain 



N 



N 



N 



- Xi) {U^+\xi)) Wi = - Xi) {Ui,{xi)) Wi-2AtJ2i^-Xi)Ui,{xi)i 



i=0 



i=0 



i=0 



N 



i=0 



dx 



(1 - Xi)Wi 



dx 



For stability we need that a certain discrete energy of be bounded in time: 



JV 



N 



- Xi) {U^+\xi)) Wi<J2i^-Xi){Ui,{xi)) Wi 



i=0 



i=0 



(174) 



which means that 



j=0 



dx 



N 



(1 - Xi)wi - 2At^(l - Xi)U^{xi)i 



dU, 



N 



i=0 



dx 



< 0. 



With the exactness of the Legendre-Gauss quadrature rule for polynomials of degree lower than 
2 A'' + 1, we have that 



N 

E 



dV^ 



dx 



„J (l-x,)». = £(^)\l-»)dx, 



and, with an additional integration by parts 



JV 



i=0 



dx 



{Ui(x)) ix. 



The stability condition obtained from energy analysis translates into an upper bound for the time- 
step, which can be seen as an accurate estimate of the CFL restriction on the time-step: 



At < 



j\{uUx))'dx 



(175) 



4.4.1 Strong stability preserving methods 

The above fully-discrete analysis must in principle be performed for every time-marching scheme. 
Therefore, it is very convenient to have a way of extending results from first-order Euler method to 
higher-order methods. Strong stability preserving Runge-Kutta and multi-step methods preserve 
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such kind of stability properties, including the case of non-linear stability analysis. A general 
review on the subject has been done by Shu [201] . and we list some results hereafter. 
If we consider the general time ODE: 

^^L^Um, (176) 

arising from the spatial discretization of the PDE (jlisp . we suppose that, after discretization in 
time using the first-order forward Euler scheme, the strong stability requirement ||f/^^'^|| < ||C^jvll 
gives a CFL restriction of the type (|175p 

At < AtpE- (177) 
We can then write an s-stage Runge-Kutta method in the form 

tAO) _ rrJ 
'-^ N — 



k=0 

and see that, as long as ai^k ^ and Pi^k ^ 0, all the intermediate stages are simply convex 
combinations of forward Euler operators. If this method is strongly stable for LjVi under the 
condition (|177p . then the intermediate stages can be bounded and the Runge-Kutta scheme is 
stable under the CFL condition 

At<cAtFE, c = min^. (178) 

In the same manner, one can devise strong stability preserving explicit multi-step methods of 
the form 

i=l 

which can also be cast into convex combinations of forward Euler steps and therefore, these multi- 
step methods are also stable, provided that 

At<cAtFE, c = min^. (179) 

i Pi 

Examples of useful coefficients for Runge-Kutta and multi-step strong stability preserving methods 
can be found in j2011 1117) . Optimal such methods are those for which the CFL coefficient c is as 
large as possible. 



4.5 Going further: High-order time schemes 

When using spectral methods in time-dependent problems, it is sometimes frustrating to have 
so accurate numerical techniques for the evaluation of spatial derivatives, and the integration of 
elliptic PDEs, whereas the time derivatives, and hyperbolic PDEs, do not benefit from the spectral 
convergence. Some tentative studies are being undertaken in order to represent also the time 
interval by spectral methods jll3| . In this spherically symmetric study of the wave equation in 
Minkowski space-time, Hennig and Ansorg have applied spectral methods to both spatial and 
time coordinates. Moreover, they have used a conformal compactification of Minkowski space- 
time, making the wave equation singular at null infinity. They have obtained nicely accurate and 
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spectrally convergent solution, even to a non-linear wave equation. If these techniques can be 
applied in general three-dimensional simulations, it would really be a great improvement. 

Nevertheless, there are other, also more sophisticated and accurate time-integration techniques 
that are currently investigated for several stiff PDEs [124] . among which Korteweg-de Vries and 
nonlinear Schrodinger equations |129j . Many such PDEs share the properties of being stiff (very 
different time-scales/ characteristic frequencies) and combining low-order non- linear terms with 
higher-order linear terms. Einstein evolution equations can also be written in such a way |37j . Let 
us consider a PDE 

— = Lu + Mu, (180) 

with the notations of Section 14.1.11 and N being a nonlinear spatial operator. Following the same 
notations and within spectral approximation, one recovers 

Fill 

—^ = LnUn+MnUn. (181) 
ot 

We detail hereafter five methods to solve this type of ODEs (see also [124] ): 

• Implicit-explicit techniques use some explicit multi-step scheme to advance the nonlinear 
part A/at, and an implicit one for the linear one. 

• Split-step are effective when the equation splits into two equations which can be directly 
integrated (see [129] for examples with the nonlinear Schrodinger and Korteweg-de Vries 
equations). 

• Integrating factor is a change of variable that allows for the exact solution of the linear 
part 

Vn = e-'^^'UN, (182) 



and to use an explicit multi-step method for the integration of the new nonlinear part 

Watc^^Vat. (183) 



dt 

• Sliders can be seen as an extension of the implicit-explicit method described above. In 
addition to splitting to problem into a linear and nonlinear part, the linear part itself is split 
into two or three regions (in Fourier space), depending on the wavenumber. Then, different 
numerical schemes are used for different groups of wavenunibers: implicit schemes for high 
wavenumbers and explicit high-order methods for the low wavenumbers. This method is 
restricted to Fourier spectral techniques in space. 

• Exponential time-differencing have been known for some time in computational electro- 
dynamics. These methods are similar to the integrating factor technique, but one considers 
the exact equation over one time-step 

= e^"^'U^ + e^"^* / e-^"WAr(i7Ar(iVAt + t), NAt + T)dT. (184) 
Jo 

Various orders for these schemes come from the approximation order of the integral. For ex- 
ample Kassam and Trefethen [124] consider a fourth-order Runge-Kutta type approximation 
to this integral, where the difficulty comes from the accurate computation of functions which 
suffer from cancellation errors. 
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5 Stationary Computations and Initial Data 



5.1 Introduction 

In this section, we restrict ourselves to problems where time does not appear explicitly. This 
is especially the case for systems which are stationary, like neutron stars in rotation or binary 
systems in circular orbits. The computation of initial data also falls into this class, given that it 
consists in finding a particular solution of Einstein equations at a given time only. Indeed, when 
using the standard 3+1 decomposition of spacetime, the initial data that are passed to the evolution 
equations cannot be totally arbitrary and must satisfy a set of equations called Einstein's constraint 
equations. For more details on the initial data problem we refer to the review by G.B. Cook [50] . 
So, in treating the problems considered here, one can forget about the issues specific to time 
presented in Sec. [H 

It must be said that spectral methods are not the only technique that has been successfully 
used to generate stationary spacetimes. The papers [Ml I222[ \6T\ 1147] give some examples of this 
fact, especially in the case of binary systems, for neutron stars or black holes. More references can 
be found in [60] . 

5.2 Single compact stars 

The computation of the structure of stationary compact stars dates back to 1939 with the famous 
solution of Tolman-Oppenheimer-Volkoff. During the last years, the need for accurate models has 
been more pressing especially with the coming online of the gravitational wave detectors which 
could help to probe the interior of such compact stars. Isolated stars in rotation are essentially 
axisymmetric but some physical effects can induce a symmetry breaking that could lead to the 
emission of gravitational waves. In the following, we shall review some computations that aim at 
including some of those effects, like the spontaneous symmetry breaking, the inclusion of magnetic 
field, the effect of exotic dense matter, mainly with strange quarks or the influence of an interior 
composed of two different superfluids. 

5.2.1 Formalisms 

The first computation of models of relativistic rotating stars in general relativity, by means of 
spectral methods, is presented in 41 . The equations are solved in spherical coordinates (see 
Section 13. 2p . Doing so, the fields only depend on the azimuthal angle 9 and on the radius r. 
The fields are expanded onto spherical harmonics with respect to the angle and onto Chebyshev 
polynomials with respect to r. The use of spherical harmonics gives a natural way of dealing 
with the coordinate singularity on the z-axis. In [41] the whole space is divided into two spherical 
domains, the outer one extending up to infinity by making use of the compactification in 1 /r seen 
in Section[3TTT2j With such setting, Einstein equations reduce to a set of four elliptic equations with 
sources extending up to infinity that are solved using a version of the algorithm based on matching 
with the homogeneous solutions (presented in Section [2.6.4^ . for each spherical harmonics. The 
system is complete once a description of the matter is given. The simplest choice is to consider a 
polytropic fiuid, with or without magnetic field. The system is solved by iteration. 

In the paper [41] , a particular emphasis is put on the various methods to measure the accuracy 
of the code. For non-rotating stars, the error is found to decrease exponentially, as the number of 
coefficients increases (see Figures 5 and 6 of [H]). However, for fast-rotating configurations, the 
error only decays as a power-law (see Figure 7 of [H]). This comes from the fact that quantities 
like the energy density are no longer C°° across the star's surface. Nevertheless, the results are in 
good agreement (to the level of a 0.1%) with those obtained by other numerical methods, as can 
be seen in [161] . 
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Spectral convergence can be recovered by using surface-adapted coordinates as first done in [36] . 
A regular mapping of the numerical coordinates to the physical ones is introduced, so that the 
surface of the star lies at the boundary between two domains (see Section [3. l.ip . For polytropes 
with 7 < 2, this is sufficient to recover spectral convergence (see Figures 5 and 6 of ^SSj). However, 
for 7 > 2, some quantities are still diverging at the surface but the convergence can be made closer 
and closer to the spectral one by analytically regularizing the density (see Section IV of [38]). 
Doing so, the error decreases as a power-law but the decrease can be made arbitrary fast at the 
cost of increasing the number of operations and so the computational time. 

Up to 2006, the neutron stars were computed using quasi-isotropic coordinates. However, in 
order to use those configurations as initial data for evolutionary codes, it may be useful to allow 
for other choices. Among the possible gauges, the Dirac one is one of the most promising [37] . 
In 1 134] models of rotating neutron stars, in the Dirac gauge are computed, for both polytropic and 
realistic equations of state. Contrary to the quasi-isotropic coordinates, the use of this new gauge 
implies to solve one tensor-like Poisson equation. Configurations obtained with the two different 
formalisms are shown to be in very good agreement. 

5.2.2 Rotating neutron star models 

Even before adapted mappings were available, interesting results could be obtained. In two pa- 
pers |1841 1185] , models of rotating neutron stars with various equations of state have been com- 
puted. Among the most surprising findings, let us mention the existence of supra-massive stars. 
Those stars do not connect to the non-rotating limit. Indeed, their high mass can only be supported 
by the presence of a centrifugal force. One of the remarkable features of such stars is the fact that 
they actually spin-up when they lose angular momentum, contrary to what is observed for normal 
stars. This effect can also be seen for neutron stars containing hyperons and thus a softer equation 
of state [238]. Let us mention that, in this case, the stability analysis of the configurations required 
the great precision that spectral methods with adapted coordinates could provide. 

It is known that isolated pulsars spin down due to magnetic braking. As the rotational frequency 
decreases, it is possible that the star will encounter a transition from one state of matter to another. 
Stationary rotating models have been used to determine the properties of such transitions |2311 
12321 1233] . A puzzling result is that the amount of energy released in a first order phase transition 
does note depend on the orbital velocity of the star and is the same as for non- rotating ones. This 
is shown to be the case for both weak [2321 and strong [233] first order transitions. 

5.2.3 Spontaneous symmetry breaking 

It is known that stars can undergo a spontaneous symmetry breaking when rotating fast enough. 
When such a phenomenon occurs, triaxial configurations are formed that are potential emitters 
of gravitational waves. The departure from axisymmetry is studied in two papers by the Meudon 
group [34] 135] . The idea of the method is to start from an axisymmetric neutron star configuration 
and to follow the growth or decay of triaxial instabilities. Well-established results in the Newto- 
nian regime are recovered and this work presents the first results in general relativity, for various 
equations of states. For a few of them, the frequency at which symmetry-breaking occurs lies in 
the frequency band of the LIGO and Virgo detectors. 

In 2002, this work has been extended [89] by making use of surface-fitting coordinates. This 
enabled the authors to obtain results in the incompressible case by properly dealing with dis- 
continuities lying at the surface of the star. Classical results in the incompressible case are thus 
recovered and it is found that the inclusion of relativity has only a moderate effect. Indeed the 
critical ratio between the kinetic energy and the absolute gravitational one T/ |W^| at which the 
triaxial instability occurs is only 30% larger for relativistic stars, with respect to their classical 
counterparts. 
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If relativistic effects only slightly stabilize the stars, the same is not true for differential rotation. 
Indeed, in [183] . the authors study various rotation profiles and equations of state using the same 
technique as in [Ml [55] to determine the onset of instability. It appears that the critical value of 
T I |W^| can be almost twice as high as for uniformly rotating stars. 

5.2.4 Configurations witii magnetic field 

Even if magnetic fields are strong in neutron stars, the structure of the objects is not affected until 
it reaches huge values, of the order of at least 10^°T. In [35], models of rapidly rotating stars 
with poloidal fields are constructed, for various equations of state. The magnetic effects are taken 
into account consistently by solving the appropriate Maxwell equations, also by means of spectral 
methods. The maximum mass of highly magnetized neutrons stars is found to be higher from 13 to 
29 % than for the non-magnetized stars. The magnetic field induces an additional pressure which 
can help to support more massive stars, thus explaining this increase. 

The presence of a magnetic field can also lead to a deformation of the neutron star. Such defor- 
mation could lead to the formation of a triaxial configuration, which would then emit gravitational 
wave. In [36] the emitted signal is computed. Typically the system radiates at two frequencies: 
and 20 where O is the angular velocity of the star. 

In a more recent work by the Meudon group [160] . magnetized configurations have been com- 
puted using coordinates matched to the surface of the star, thus making the computation much 
more accurate. Gyromagnetic ratios of rapidly rotating neutron stars of various equations of state 
are obtained. The limit of a ratio g = 2, value for a charged black hole, is never reached. 

5.2.5 Strange stars 

It is possible that the fundamental state of nuclear matter is not the ordinary matter but rather 
a plasma of deconfined quarks u, d and s, called strange matter. If this is the case, neutron stars 
would rather be strange stars. The main difference between those two types of compact stars 
is that strange ones are somewhat smaller and thus more compact. In particular, they would 
support higher rotation rates. There is a strong density jump at the surface of a strange star and 
surface-fitting coordinates are required in order to deal with it. 

Fast rotating models of strange stars are computed in fl031 188]. Due to higher compactness, it 
is found that strange stars can rotate significantly faster than their neutron star counterparts. The 
attained T / \W\ can be twice as large. As in the neutron star case, supermassive configurations 
that spin-up with angular momentum loss are found. The influence of strange matter on the 
emission of gravitational waves is studied in [90] where viscosity effects and triaxial instabilities 
are carefully taken into account. 

It is believed that millisecond pulsars have been spun-up by accreting matter from a companion. 
However, the details of this mechanism depend on the nature of the compact object. In [237] . the 
differences between accretion onto a neutron star and onto a strange star are investigated, using 
2D stationary models computed by spectral methods. 

5.2.6 Quasi-periodic oscillations 

Quasiperiodic oscillations (QPOs) are observed in the kHz regime and are believed to be the 
signature of matter falling onto a compact object. In the standard picture, the frequency of the 
QPOs, is that of the last stable orbit around the compact object. Let us mention that the presence 
of a last stable orbit around an extended body is not an effect of relativity but can also be seen in 
the Newtonian regime, as shown in [234] . 
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The precise structure of the accreting object has a great influence on the QPO. In a series 
of papers |235l [521 131 HZ] , comparisons are made between observations and various compact stars 
models that could account for QPOs. 

Using a multi-domain approach, strange stars with a crust can also be computed [236] . one 
domain describing the interior of the star and another one the crust. It is shown that the presence 
of the crust could change the value of the QPO by up to 20%. 

5.2.7 More complex configurations 

In this section, objects in more exotic configurations are presented. This is an illustration of 
both the complexity of the neutron stars physics and the ability of spectral methods to deal with 
complicated systems. 

The observation of glitches in isolated pulsars is consistent with the presence of a superfluid 
interior. The simplest model considers two fluids, one composed of neutrons and the other one of 
protons and electrons, both components being superfluids. However, those two components could 
have different velocities, in particular different rotation rates. Such configurations are computed 
in [178] . A multi-domain setting is crucial to be able to follow the two different fluids because the 
components do not have the same shape. Among the various results obtained, let us mention the 
confirmation of the existence of prolate-oblate configurations. 

Neutron stars are usually assumed to be at zero-temperature. However this approximation is 
no longer true for newborn neutron stars, just after the supernova. Models of newborn neutron 
stars in uniform rotations are constructed in |106j using an extension of the code developed in [41] . 
Various hypothesis about the interior (different lepton numbers, isothermal versus isentropic) are 
considered. Sequences of fixed baryon mass and angular momentum are constructed. Such se- 
quences approximate the evolution of the proto-neutron star into a cold neutron star. The results 
have been extended to differentially rotation proto-neutron stars in [107] . 

The effect of finite temperature is also investigated in |226| . The authors found that newborn 
neutron stars are unlikely to undergo the bar mode instability but that the secular ones could take 
place and result in a significant emission of gravitational waves. Another interesting result of [226] 
is the existence of toroidal-like configurations, which appear for a broad range of parameters and 
before the mass-shedding limit. In such cases, the surface of the star is highly deformed and 
surface-fitting coordinates are required. 

Axisymmetric rotating neutron stars have also been computed by a code developed by M. An- 
sorg and collaborators [21 [10] . This code is based on Lewis-Papapetrou coordinates (p,^) , which 
are closely related to the usual cylindrical coordinates. Typically, space is decomposed into two 
domains: one for the interior of the star and another to the exterior which extends up to spatial 
infinity. Compactification of space and surface-fitting mappings are used. Both variables are ex- 
panded on Chebyshev polynomials. Instead of solving the equations harmonics by harmonics and 
iterate, as is done by the Meudon group, the equations are written with a collocation method (see 
Section I2.5.3|) and solved as one single system. The price to pay is that the size of the system is 
somewhat larger (i.e. in m^, m being the total number of coefficients for each coordinates). The 
system is solved by means of the Newton-Raphson's method. At each step, the linear system is 
solved using iterative techniques with preconditioning. With this setting, impressive accuracy is 
reached. 

The coordinates used in |9l [10] are more general than the ones used by the Meudon group, 
especially with respect to their surface-fitting capabilities. They can account for more complicated 
configurations and, in particular, highly distorted matter distribution can be obtained. This is 
shown in [121 111] , where relativistic axisymmetric toroidal configurations of matter, known as the 
Dyson rings, are computed. Such rings are obtained up to the mass-shedding limit. Transition to 
the limit of an extreme Kerr black hole is also discussed. 
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5.3 Single black holes 

Compared to the compact star case, single black holes have not been very much studied. This is 
probably because the structure of a stationary black hole is somewhat simpler than the one of a 
compact star. However, as it will be seen, there are still properties that must be investigated. 

Spacetimes containing a single black hole constitute a good benchmark for numerical methods, 
a lot of results being known analytically. In [125] , the authors have implemented a spectral solver 
and applied it to various test problems. The solver itself is two-dimensional and thus applicable 
only to axisymmetric systems. There is a single domain that consists of the whole space outside a 
sphere of given radius (i.e. the black hole). Space is compactified by using the standard variable 
l/r. The two physical variables (r,d) are mapped onto squares in and then expanded on 
Chebyshev polynomials. The equations are written using a 2-dimensional collocation method (see 
Section I2.5.3P and the resulting system is solved by an iterative algorithm (here Richardson's 
method with preconditioning). This solver is applied to solve the Einstein's constraint equations 
for three different systems: i) a single black hole ii) a single black hole with angular momentum 
iii) a black hole plus Brill waves. In all three cases, spectral convergence is achieved and accuracy 
of the order of lO""'^" is reached with 30 points in each dimension. 

A black hole is somewhat simpler than a neutron star, in the sense that there is no need for a 
description of matter (no equation of state for instance). However, in some formalisms, the presence 
of a black hole is enforced by imposing non-trivial solution on some surfaces (typically spheres). 
The basic idea is to demand that the surface is a trapped surface. Such surfaces are known to lie 
inside event horizons and so are consistent with the presence of a black hole. Discussions about such 
boundary conditions can be found in '62\- However, in non-stationary cases, the set of equations to 
be used is not easy to derive. The authors of [123] implemented various sets of boundary conditions 
to investigate their properties. Two different and independent spectral codes are used. Both codes 
are very close to those used in the case of neutron stars, one of them being based on Lorene 
library [9^ (see Section [5.2. ip and the other one has been developed by M. Ansorg and shares a 
lot a features with [9l \10\ . Such numerical tools have proved useful in clarifying the properties of 
some sets of boundary conditions that could be imposed on black hole horizons. 

The reverse problem is also important in the context of numerical relativity. In some cases one 
needs to know if a given configuration contains a trapped surface and if it can be located, at each 
time-step. Several algorithms have been proposed in the past to find the locus where the expansion 
of the outgoing light rays vanishes (thus defining the trapped surface). Even if the term is not 
used explicitly, the first application of spectral expansions to this problem is detailed in [23] . The 
various fields are expanded on a basis of symmetric trace-free tensors. The algorithm is applied 
to spacetimes containing one or two black holes. However results seem to indicate that high order 
expansions are required to locate the horizons with a sufficient precision. 

More recently, another code '135' using spectral method has been used to locate apparent 
horizons. It is based on the Lorene library with its standard setting i.e. a multi-domain de- 
composition of space and spherical coordinates (see Section 15.2.11 for more details). The horizon 
finder has been successfully tested on known configurations, like Kerr-Schild black holes. The use 
of spectral methods makes it both fast and accurate. Even if the code is using only one set of 
spherical coordinates (hence its presentation in this section), it can be applied to situations with 
more than one black hole, like the well-known Brill-Lindquist data [50] . 

5.4 Rings around black holes 

The problem of uniformly rotating rings surrounding a black hole can be viewed as an intermediate 
step between one body, axisymmetric configurations and the two body problem. Indeed, even if 
one has to deal with two components, the problem is still axisymmetric. In 13 , configurations of 
a black hole surrounded by a uniformly rotating ring of matter are computed in general relativity. 
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The matter is assumed to be a perfect fluid. To solve the equations, space is divided into five 
computational domains. One of them describes the ring itself, another one the region around the 
black hole and another extends up to infinity. One of the difficulties is that the surface of the ring is 
not know a priori and so the domains must be dynamically adapted to its surface. Cylindrical-type 
coordinates are used and, in each domain, are mapped onto squares of numerical coordinates. The 
actual mappings depend on the domain and can be found in Section IV of [13] . 

Numerical coordinates are expanded onto Chebyshev polynomials. The system to be solved is 
obtained by writing Einstein equations in the collocation space including regularity conditions on 
the axis and appropriate boundary conditions on both the horizon of the black hole and at spatial 
infinity. As in [9l[l0], the system is solved iteratively, using Newton-Raphson's method. 

Both Newtonian and relativistic configurations are computed. The ratio between the mass of 
the black hole and the mass of the ring is varied from (no black hole) up to 144. The inner 
mass shedding of the ring can be obtained. One of the most interesting results is the existence of 
configurations for which the ratio Jc/A/^ of the black hole angular momentum and the square of 
its mass exceeds one, contrary to what can be achieved for an isolated black hole. 

5.5 Binary compact stars 
5.5.1 Formalism 

Systems consisting of two compact objects are known to emit gravitational waves. Due to this 
emission, no closed orbits can exist and the objects follow a spiral-like trajectory. It implies that 
such systems have no symmetries that can be taken into account and full time-evolutions should 
be made. However, when the objects are relatively far apart, the emission of gravitational waves is 
small and the inspiral slow. In this regime, one can hope to approximate the real trajectory with 
a sequence of closed orbits. Moreover, the emission of gravitational waves is known to efficiently 
circularize eccentric orbits so that only circular orbits are usually considered. 

So, a lot of efforts have been devoted to the computation of circular orbits in general relativity. 
This can be done by demanding that the system admit an helical Killing vector dt -f ^^^p , fl being 
the orbital angular velocity of the system. Roughly speaking, this means that advancing in time 
is equivalent to turning the system around its axis. Working in the corotating frame, one is left 
with a time-independent problem. 

Additional approximations must be made in order to avoid any diverging quantities. Indeed, 
when using the helical symmetry, the system has an infinite lifetime and can fill the whole space 
with gravitational waves, thus causing quantities like the total mass to be infinite. To deal with 
that, various techniques can be used, the simplest one consisting in preventing the appearance of 
any gravitational waves. This is usually done by demanding that the spatial metric be conformally 
flat. This is not a choice of coordinates but a true approximation that has a priori no reason to be 
verified. Indeed, even for a single rotating black hole, one can not find coordinates in which the 
spatial 3- metric is conformally flat, so that we do not expect it to be the case for binary systems. 
However comparisons with post-Newtonian results or non-conformally flat results tend to indicate 
that this approximation is relatively good. 

Under these assumptions, Einstein equations reduce to a set of flve elliptic equations for the 
lapse, the conformal factor and the shift vector. Those equations encompass both the Hamiltonian 
and momentum constraint equations and the trace of the evolution equations. To close the system, 
one must provide a description of the matter. It is commonly admitted that the fluid is irrotational, 
the viscosity in neutron stars being too small to synchronize the motion of the fluid with the orbital 
motion. It follows that the motion of the fluid is described by an additional elliptic equation for 
the potential of the flow. The matter terms entering the equations via the stress-energy tensor can 
then be computed, once an equation of state is given. An evolutionary sequence can be obtained 
by varying the separation between both stars. 
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5.5.2 Numerical procedure 

Up to now, only the Meudon group has solved those equations by means of spectral methods in 
the case of two neutron stars. Two sets of domains are used, one centered on each star. Each 
set consists of spherical-like domains that match the surface of the star and extend up to infinity. 
Functions are expanded onto spherical harmonics with respect to the angles [9, ip) and Chebyshev 
polynomials with respect to the radial coordinates. Each Poisson equation AA^ = Sn is split into 
two parts ATVi = Sni and AN2 = Sn2, such that Sn = Sni + and N = Ni + N2. The 
splitting is of course not unique and only requires that Sn- is mainly centered around the star i so 
that it is well described by spherical coordinates around it. The equation labeled i is then solved 
using the domains centered on the appropriate star. The splittings used for the various equations 
can be found explicitly in Section IV-C of [lOOj . 

The elliptic equations are solved using the standard approach by the Meudon group found 
in [109] . For each spherical harmonic, the equation is solved using a Tau-method and the matching 
between the various domains is made using the homogeneous solutions method (see Section [2. 6. 4[) . 
The whole system of equations is solved by iteration and most of the computational time is spent 
when quantities are passed from on set of domains to the other one by means of a spectral sum- 
mation (this requires operations, N being the number of collocation points in one dimension). 
A complete and precise description of the overall procedure can be found in [lOOj . 

5.5.3 Binary neutron stars 

The first sequence of irrotational neutron star binaries computed by spectral means is shown 
in [40]. Both stars are assumed to be polytropes with an index j — 2. The results are in good 
agreement with those obtained, simultaneously, by other groups with other numerical techniques 
(see for instance [2411222] ). One of the important points that has been clarified by [40] concerns the 
evolution of the central density of the stars. Indeed, at the end of the nineties, there was a claim 
that both stars could individually collapse to black holes before coalescence, due to the increase 
of central density as the two objects spiral towards each other. Should that have been true, this 
would have had a great impact on the emitted gravitational wave signal. However it turned out 
that this was coming from a mistake in the computation of one of the matter term. The correct 
behavior, confirmed by various groups and in particular by [40] , is a decrease in the central density 
as the stars get closer and closer (see Figure I of [40]). 

If the first sequence computed by spectral methods is shown in [4^ , the complete description and 
validation of the method are given in [100] . Convergence of the results with respect to the number 
of collocation points is exhibited. Global quantities like the total energy or angular momentum 
are plotted as a function of the separation and show remarkable agreement with results coming 
from analytical methods (see Figures 8 to 15 of [100] ). Relativistic configurations are also shown 
to converge to the Newtonian counterparts when the compactness of the stars is small (Figures 16 

to 20 of [Ton] ). 

Newtonian configurations of compact stars with various equations of state are computed for 
both equal masses [216] and various mass ratios [213] . One of the main results of the computations 
concerns the nature of the end point of the sequence. For equal masses, the sequence ends at 
contact for synchronized binaries and at mass shedding for irrotational configurations. This is to 
be contrasted with the non-equal mass case where the sequence always ends at the mass shedding 
limit of the smallest object. 

Properties of the sequences in the relativistic regime are discussed in [2141 1215] . In [214] 
sequences with j — 2 are computed, for various compactness and mass ratios, for both synchronized 
and irrotational binaries. The nature of the end point of the sequences is discussed and similar 
behavior to the Newtonian regime is observed. The existence of a configuration of maximum 
binding energy is also discussed. Such existence could have observational implications because 
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it could be an indication of the onset of a dynamical instability. Sequences of polytropes with 
various indexes ranging from 1.8 to 2.2 are discussed in |215j . In particular, the authors are lead 
to conjecture that, if a configuration of maximum binding energy is observed in Newtonian regime, 
it is also observed in conformal relativity for the same set of parameters. 

In [76] the authors derive, from the sequences computed in [214] . a method to constrain the 
compactness of the stars from the observations. Indeed, from results |214j . one can easily determine 
the energy emitted in gravitational waves per interval of frequency (i.e. the power-spectrum of the 
signal). For large separation, that is for small frequencies, the curves follow the Newtonian one. 
However, there is a break frequency at the higher end (see Figure 2 of [76]). The location of this 
frequency depends mainly on the compactness of the stars. More precisely, the more compact the 
stars are, the higher the break frequency is. Should such frequency be observed by the gravitational 
wave detectors, this could help to put constraints on the compactness of the neutron stars and 
thus on the equation of state of such objects. 

5.5.4 Extensions 

The framework of [100] is applied to more realistic neutron stars in [26] . In this work, the equations 
of state are more realistic than simple polytropes. Indeed, three different equations are considered 
for the interior, all based on different microscopic models. The crust is also modeled. For all the 
models, the end point of the evolution seems to be given by the mass shedding limit. However, 
the frequency at which the shedding occurs depends rather strongly on the EOS. The results are 
in good agreement with post-Newtonian ones, until hydrodynamic effects begin to be dominant. 
This occurs at frequencies in the range 500-1000 Hz, depending on the EOS. 

Sequences of binary strange stars have also been computed [133j . Contrary to the neutron star 
case, the matter density does not vanish at the surface of the stars and one really needs to use 
surface-fitting domains to avoid any Gibbs phenomenon that would spoil the convergence of the 
overall procedure. Sequences are computed for both synchronized and irrotational binaries and a 
configuration of maximum binding energy is attained in both cases. This is not surprising: the 
strange stars are more compact than the neutron stars and are less likely to be tidally destroyed 
before reaching the extremum of energy, making it easier to attain dynamical instability. More 
detailed results on both neutron star and strange star binaries are discussed in the follow-up 
papers [871I9T] . 

All the works presented above are done in the conformal flatness approximation. As already 
stated in Section 15.5.11 this is only an approximation and one expects that the true conformal 
3-metric will depart from flatness. However, in order to maintain asymptotic flatness of space- 
time, one needs to get rid of the gravitational wave content. One such waveless approximation 
is presented in |196| and implemented in [223 . Two independent codes are used, one of them 
being an extension of the work described in [100] . The number of equations to be solved is then 
greater than in conformal flatness (one has to solve for the conformal metric), but the algorithms 
are essentially the same. It turned out that the deviation from conformal flatness is rather small. 
The new configurations are slightly further from post-Newtonian results than the conformally flat 
ones, which is rather counter-intuitive and might be linked to a difference in the definition of the 
waveless approximations. 

5.6 Binary black hole systems 
5.6.1 Digging the holes 

If the computation of binary black holes in circular orbits has a lot of common features with the 
neutron star case, there are also some differences that need to be addressed. In at least one aspect, 
black holes are much simpler objects because they are solution of Einstein equations without 
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matter. So the whole issue of investigating various equations of state is irrelevant and there is no 
need to solve any equation for the matter. However, there is a price to pay and one must find a 
way to impose the presence of holes in the spacetime. Two main ideas have been proposed. 

In the puncture method, the full spacetime contains three asymptotically flat regions. One is 
located at r = cx) and the two others at two points Ml an M2 which are called the punctures. 
The presence of the flat regions near the punctures is enforced by demanding that some quantities, 
like the conformal factor, diverge at the those points (typically like l/r). The discontinuities are 
taken out analytically and the equations are solved numerically for the regular parts, in the whole 
space. This idea dates back to the work of Brill and Lindquist '50], at least in the case of black 
holes initially at rest. The puncture approach has been successfully applied to the computation of 
quasi-circular orbits by means of spectral methods in [8]. 

The apparent horizon method relies on initial works by Misner [151 ' and Lindquist [140' . In this 
case, the space has only two asymptotically flat regions. One can show that this is equivalent to 
solving Einstein's equations outside two spheres on which boundary conditions must be imposed. 
The boundary conditions are based on the concept of trapped surface and apparent horizons. The 
physical state of the black holes are precisely encoded in the boundary conditions. 

5.6.2 First configurations 

The first configurations of binary black holes computed by means of spectral methods can be 
found in |110] . The formalism and various hypothesis are given in the companion paper [98]. The 
assumptions are very similar to those used for binary neutron stars (see Section I5.5.1|) . Helical 
symmetry is enforced and conformal flatness assumed. The holes are described by the apparent 
horizon technique. However, the boundary conditions used have been shown to be only approxi- 
mately valid, up to a rather good accuracy. This effect is discussed in the original paper [110] and 
further explored by Cook in [62] . The numerical techniques are very similar to the ones employed 
for binary neutron star configurations (see Section [5. 5. 2p . Two sets of spherical domains are used, 
one for each black hole. Boundary conditions are imposed on the surface between the nucleus and 
the first shell. Both sets extend up to infinity using a compactification in 1/r. 

For the first time, a good agreement was found between numerical results and post-Newtonian 
ones. A detailed comparison can be found in |67| . In particular, the location of the minimum 
of energy is shown to coincide at the level of a few percent. This improvement with respect to 
previous numerical works is mainly due to a difference in the physical hypothesis (i.e. the use 
of helical symmetry) . One important product of [110] is the use of a new criterion to determine 
the appropriate value of the orbital angular velocity J7. Indeed, for neutron stars, this is done by 
demanding that the fluid of both stars be in equilibrium [100) . This, of course, is not applicable for 
black holes. Instead, in [98l 1110) . it is proposed to find by demanding that the ADM mass and 
the Komar-like mass coincide. Without going into to much details, this amounts to demanding 
that, far from the binary and at first order in 1/r, the metric behave like the Schwarzschild one. It 
is shown in [98] that it can be linked to a relativistic virial theorem. Since then, it has been shown 
that this criterion could also be used for neutron stars [215] and that it was equivalent to the use 
of a variational principle called the effective potential method [611 1221 11731 159] , where the binding 
energy is minimized with respect to fi. 

5.6.3 Further investigations 

More recently, two other spectral codes have been developed in the context of binary black holes 
and successfully applied to address some of the issues raised by the work of [HHl HlOj . 

One of those codes is due to the Caltech/Cornell group by H. Pfeiffer and collaborators and 
is described extensively in [1721 1168) . The code is multi-domain and two main types of domains 
are used i) square domains where each Cartesian-like coordinate is expanded onto Chebyshev 
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polynomials and ii) spherical domains where spherical harmonics are used for the angles (0, (p) and 
Chebyshev polynomials for the radial coordinate. Space can be compactified by a standard use 
of the variable l/r. The two types of domains can be set up in various manners to accommodate 
the desired geometry. When using both square and spherical domains, there are regions of space 
that belong to more than one domain. This is to be contrasted with work by the Meudon group 
where domains are only touching but not overlapping. The algorithm of [172] solves differential 
equations by using a multi-dimensional collocation method. The size of the resulting system is 
roughly equal to the number of collocation points. It is then solved iteratively via a Newton- 
Raphson algorithm with line search. At each step of the Newton-Raphson method, the linear 
system is solved by means of an iterative scheme (typically GMRES). This inner iterative solver 
requires careful preconditioning to work properly. Various tests are passed by the code in [172] . 
where elliptic equations and systems are solved in cither spherical or bispherical topologies. In the 
cases presented the error decays spectrally. 

In 1 170] the code is used to investigate different ways of solving the constraint equations. 
Three different decompositions are used: the conformal TT one, the physical TT one and the 
thin-sandwich decomposition. When solving for the constraint equations only, one also needs to 
precise some freely specifiable variables, which describe the physical state of the system. In [170] . 
those specifiable variables are fixed using a superposition of two Kerr-Schild black holes. The net 
result of [170] is that global quantities, like the total energy, are very sensitive to the choice of 
decomposition. The variation of total energy can be as large as 5%, which is the order of the energy 
released by gravitational waves. It is also shown that the choice of extrinsic curvature tensor is 
more crucial than the one of conformal metric, in accordance with an underlying result of [110] . 
Let us precise that the equations derived form the helical Killing vector approach in [98j 1110] 
are equivalent to the ones obtained by making use of the thin-sandwich decomposition of the 
constraints. The freely specifiable variables are obtained by both the imposition of the helical 
Killing symmetry and by solving an additional equation for the lapse function (resulting in the 
so-called extended thin-sandwich formalism). 

In ^63j the boundary conditions based on the apparent horizon formalism [62] are implemented 
and tested numerically in the case of one or two black holes. In the latter case, the main difference 
with 110 , lies in the use of more elaborate and better boundary conditions on the horizons of the 
holes. By allowing for non-vanishing lapse on the horizons, the authors of 63 solve the constraint 
equations exactly. This is to be contrasted with [110] . where the momentum constraint equation 
was only solved up to some small correction. Anyway, both results show a rather good agreement. 
This is not surprising because the correction used by the Meudon group was known to be small (see 
Figures 10 and 11 of [110] ). More results are presented in 59 , for both corotating and irrotational 
black holes. An important result of [591 is the comparison of the two criteria for determining the 
orbital angular velocity Q. They indeed show that the the effective potential method first introduced 
in [61] and the method based on the virial theorem proposed in [98] are in very good agreement. 

By slightly extending the boundary conditions used in [59], the authors of [169^ proposed to 
reduce the eccentricity of the binary black hole configurations. This is done by giving the holes 
a small radial velocity by modifying the boundary condition on the shift vector. The code and 
other equations are the same as in [59] . Time evolution of the obtained initial data shows indeed 
that this technique can reduce the eccentricity of the binary. However, the effect on the emitted 
gravitational wave is small and probably unimportant. 

Another application of the Caltech/Cornell solver can be found in |143] . where the emphasis 
is put on nearly maximum spinning black holes. Initial data are constructed for both single and 
binary black holes. Three families of initial data are investigated. Using a formalism based on the 
Kerr-Schild spacetime, the authors are able to reach spins as large as a = 0.9998. Such nearly 
maximum spinning black holes may be relevant from the astrophysical point of view. Evolutions 
of those data are also discussed. 
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The other spectral code used to compute configuration of binary black holes is due to M. An- 
sorg [5]. It shares a lot of features with the code developed by the same author in the context 
of rotating stars [51 HO] already discussed in Section 15.2.71 Space is decomposed in two domains. 
One of them lies just outside the horizons of the holes and bispherical-like coordinates are used. 
The other domain extends up to infinity and an appropriate mapping is used in order to avoid 
the singularity of the bispherical coordinates at spatial infinity (see Section IV of [6]). The angle 
of the bispherical coordinates (i.e. the angle around the x-axis joining the two holes) is expanded 
onto Fourier series and the two other coordinates onto Chebyshev polynomials. Like in I172| . 
the partial differential equations are solved using a collocation method and the resulting system is 
solved by Newton-Raphson's method. Once again the linear system coming from the Jacobian is 
solved by an iterative scheme with preconditioning. The code is used to compute essentially the 
same configuration as those shown in [59j . An interesting point of [6] is the detailed investigation 
of convergence of the results when increasing the resolution. As can bee seen in Figure 4 of [Hj, 
the error starts by decreasing exponentially but, for high number of points, it seems that the error 
only follows a power-law. This is an indication that some non-C°° fields must be present. It is 
conjectured in [6^ that this comes from logarithm terms that can not be dealt with properly with 
a compactification in 1/r. The same kind of effect is investigated in some details in |109j . where 
some criteria for the appearance of such terms are discussed. 

A code very similar to the one used in |6J has also been used to compute spacetimes with black 
holes using the puncture approach [8] . Given that the black holes are no longer described by their 
horizons, one do not need to impose inner boundary conditions. The absence of this requirement 
enabled the author of |8] to use a single domain to describe the whole space, from the puncture 
up to infinity. The other features of the spectral solver are the same as in [6 . This method has 
been successfully applied to the computation of binary black hole configurations in the puncture 
framework. The authors have, in particular, investigated high mass ratios between the bodies and 
compared their results with the ones given in the test-mass limit around a Schwarzschild black 
hole. The discrepancy is found to be of the order of 50% for the total energy. It is believed that 
this comes from the fact that the mass of each puncture cannot be directly related to the local 
black hole mass (see discussion in Section VII of [S]). 

Let us finally mention that the algorithms developed by M. Ansorg in [9l (TOj [H [6] have all been 
unified in [7j to accommodate any type of binaries. Various domain decompositions are exhibited 
that can be used to represent neutron stars, excised black holes or puncture black holes, with 
compactification of space. The algorithms are shown to be applicable to limiting cases like large 
mass ratios. 

5.7 Black hole- neutron star binaries 

Until recently the binaries consisting of a neutron star and a black hole received fewer attention 
than the other types of systems. It was believed, and this was partly true, that this case could 
easily be handled once the cases of binary neutron stars and binary black holes were understood. 
However, such binaries are of evident observational interest and could be the most promising source 
of gravitational waves for the ground-based detectors [28] . 

The first application of spectral methods to the black hole-neutron star binaries can be found 
in [209] . The main approximation is to consider that the black hole is not influenced by the neutron 
star. Technically, this means that the Einstein equations are split into two parts (i.e. like for binary 
neutron stars I5.5.2P . However the part of the fields associated to the black hole are fixed to their 
analytical value. As the fields are not solved for the black hole part, the results should depend 
on the actual splitting, the equations being non-linear. The part of the fields associated with the 
neutron star are solved using the standard setting for the Meudon group. Of course, this whole 
procedure is only valid if the black hole is much more massive than the neutron star and this is 
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why |209) is limited to mass ratios of 10. In this particular case, it is shown that the results depend 
to the level of a few percent on the choice of splitting, which is a measure of the reached accuracy. 
It is also shown that the state of rotation of the star (i.e. synchronized or irrotational) has little 
influence on the results. 

In [210] the equations of the extended thin-sandwich formulation are solved consistently. Like 
for the binary neutron star case, two sets of spherical coordinates are used, one centered around each 
object. The freely specifiable variables are derived from the Kerr-Schild approach. Configurations 
are obtained with a moderate mass ratio of 5. However the agreement with post-Newtonian results 
is not very good and the data seem to be rather noisy (especially the deformation of the star). 

Quasi-equilibrium configurations based on a helical Killing vector and conformal flatness have 
been obtained independently by [108] and [211) . Both codes are based on the Lorene library 
and use two sets of spherical coordinates. They differ mainly in the choice of boundary conditions 
for the black hole. However, it is shown in the erratum of 108 that the results match pretty well 
and are in very good agreement with post-Newtonian results. Mass ratios ranging from 1 to 10 
are obtained in [211] and the emitted energy spectra are estimated. The work of [211 has been 
extended in [2121 where the parameter space of the binary is extensively explored. In particular, 
the authors determine whether the end-point of the sequences is due to an instability or to the 
mass-shedding limit. It turns out that the star is more likely to reach the mass-shedding limit if it 
is less compact and if the mass ratio between the black hole and the star is important, as expected. 

More recently, the Caltech/CorncU group has applied the spectral solver of |172[ 1168] in order 
to compute black hole- neutron stars configurations [80] . Some extensions have been made to 
enable the code to deal with matter by making use of surface fitting coordinates. Thanks to the 
domain decomposition used (analogous to the one of [172[ 1168] ). the authors of [80] can reach an 
estimated accuracy of 5 • 10"^, which is better than the precision of previous works (by roughly an 
order of magnitude). Configurations with one spinning black hole and configurations with reduced 
eccentricity are also presented, in the line of |169j . 

5.8 Spacetimes with waves 

The work [171! presents a method to produce initial data configuration containing waves. Given 
a fixed background metric, it shows how to superimpose a given gravitational wave content. The 
equations are solved numerically using a multi-domain spectral code based on [1721 1168) . Space 
is covered by various spherical-like shells and is described up to infinity. When no black hole is 
present, the origin is covered by a square domain because regularity conditions at the origin, in 
spherical coordinates, are not handled by [172[ 1168] . Such setting is used to generate spacetimes 
containing i) pure quadrupolar waves ii) flat space with ingoing pulse and iii) a single black hole 
superimposed with an ingoing quadrupolar wave. 

5.9 Hyperboloidal initial data 

If the 3-1-1 decomposition is the most widely used for numerical relativity, some other possibilities 
have been proposed, with possibly better features. In particular, one can vary the foliation of 
spacetime to get hyperboloidal data. With such a setting, at infinity spacetime is foliated by light 
cones instead of spatial hypersurfaces, which makes extraction of gravitational waves in principle 
easier. 

In [81] Frauendiener is interested in generating hyperboloidal initial data sets from data in 
physical space. The technique proceeds in two steps. First a non-linear partial differential equation 
(the Yamabe equation) must be solved to determine the appropriate conformal factor tu. Then, 
the data are constructed by dividing some quantities by this w. This second step involves an 
additional difficulty: u vanishes at infinity but the ratios are finite and smooth. It is demonstrated 
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in [81] that spectral methods can deal with those two steps. Some symmetry is assumed so that 
the problem reduces to a 2-dimensional one. The first variable is periodic and expanded onto 
Fourier series whereas Chebyshev polynomials are used for the other one. The Yamabe equation 
is solved using an iterative scheme based on Richardson's iteration procedure. The construction of 
the fields, hence the division by a field vanishing at infinity, is then handled by making use of the 
non-local nature of the spectral expansion (i.e. by working in the coefficient space; see Section 4 
of [H] for more details). 
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6 Dynamical Evolution of Relativistic Systems 



The modeling of time-dependent physical systems is traditionally the ultimate goal in numerical 
simulation. Within the field of numerical relativity, the need for studies of dynamical systems 
is even more pronounced because of the seek for gravitational wave patterns. Unfortunately, as 
presented in Section 14.11 there is no efficient spectral time discretization yet and one normally 
uses finite-order time-differentiation schemes. Therefore, in order to get high temporal accuracy, 
one needs to use high-order explicit time marching schemes (e.g. fourth or sixth-order Runge- 
Kutta [49] ) . This requires quite some computational power and might explain why, except for 
gravitational collapse [§51 1157] , very few studies using spectral methods have dealt with dynamical 
situations until the Caltech/Cornell group began to use spectral methods in numerical relativity, in 
the beginning of years 2000 j l28] 1127] . This group now have a very- well developed pseudo-spectral 
collocation code "Spectral Einstein Code" (SpEC), for the solution of full three-dimensional dy- 
namical Einstein equations. 

In this section, we review the status of the numerical simulations using spectral methods in some 
fields of General Relativity and Relativistic Astrophysics. Although we may give at the beginning 
of every section a very short introduction to the context of the relevant numerical simulations, 
it is not our point to detail them since dedicated reviews exist for most of the themes presented 
here and the interested reader should consult them for physical details and comparisons with other 
numerical and analytical techniques. Among the systems which have been studied, one can find 
gravitational collapse [84] (supernova core collapse or collapse of a neutron star to a black hole), 
oscillations of relativistic stars [2071 1130] and evolution of "vacuum" spacetimes. These include 
the cases of pure gravitational waves or scalar fields, evolving in the vicinity of a black hole or as 
(self-gravitating) perturbations of Minkowski flat spacetime. Finally, we shall discuss the situation 
of compact binaries |1751 131] spectral numerical simulations. 

6.1 Single Stars 

The numerical study of the evolution of stars in General Relativity involves two parts: first one 
has to solve for the evolution of matter (relativistic hydrodynamics, see [77]), and second one 
must compute the new configuration of the gravitational field. Whereas, spectral-methods based 
codes are now able to study quite well the second part (see Section 16. 2p , the first part has not 
benefited from so many efforts from the groups using spectral methods in the past decade. One is 
facing the paradox: spectral methods have been primarily developed for the simulation of hydro- 
dynamic systems (see Section FOj) but they are not often used for relativistic hydrodynamics. This 
might be understood as a consequence of the general problem of spectral methods to deal with 
discontinuous fields and supersonic fiows: the Gibbs phenomenon (see Section [2.4.4]) . Relativistic 
flows in astrophysics are often supersonic and therefore contain shocks. Although some techniques 
have been devised to deal with them in one-dimensional studies (see e.g. |45]). there have been 
no multi-dimensional convincing work. Other problems coming from multi-dimensional relativistic 
hydrodynamics which can spoil the rapid convergence properties are the density sharp profiles 
near neutron star surfaces. These can imply a diverging or discontinuous radial derivative of the 
density, thus slowing down the convergence of the spectral series. 

6.1.1 Supernova core collapse 

The physical scenario studied here is the formation of a neutron star from the gravitational collapse 
of degenerate stellar core. This core can be thought as to be the iron core of a massive star at 
the end of its evolution (standard mechanism of type II supernova). The degeneracy pressure of 
the electrons can no longer support the gravity and the collapse occurs. When the central density 
reaches nuclear values, the strong interaction stiffens the equation of state, stopping the collapse 
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in the central region and a strong shock is generated. This mechanism has been long thought to 
be a powerful source of gravitational radiation, but recent simulations show that the efficiency is 
much lower than previously estimated [701 1196] . The first numerical study of this problem was 
the spherically symmetric approach by May and White |149j . using artificial viscosity to damp 
the spurious numerical oscillations caused by the presence of shock waves in the flow solution. 
Nowadays, state-of-the-art codes use more sophisticated High-Resolution Shock-Capturing (HRSC) 
schemes or High-Resolution Central (HRC) Schemes (for details about these techniques, see the 
review by Font [77]). The first axisymmetric fully (general) relativistic simulations of the core 
collapse scenario have been done by Shibata [192) . Shibata and Sekiguchi [196] . which then used 
HRSC schemes and a parametric equation of state. More recently, magnetohydrodynamic effects 
have been taken into account in the axisymmetric core collapse by Shibata et al. [194] . using HRC 
schemes. Three-dimensional core collapses simulations, including more realistic equation of state 
and deleptonization scheme have been performed within the CACTUS-CARPET-whisky [1481 [T7] 
framework by Ott et al. [1651 1164] . These simulations have been compared with those of the 
CoCoNuT code (see hereafter). A more detailed historical presentation can be found in the 
Living Review by Fryer and New [84] . 

The appearance of a strong hydrodynaniic shock is, in principle, a serious problem to numerical 
models using spectral methods. Nevertheless, a first preliminary study in spherical symmetry and 
in Newtonian theory of gravity has been undertaken in 1986 by Bonazzola and Marck ^43: , with 
the use of "natural" viscosity. The authors showed a mass conservation to a level better than 10"'' 
using one domain with only 33 Chebyshev polynomials. In 1993, the same authors performed the 
first three-dimensional simulation (still in Newtonian theory) of the pre-bounce phase |46j . giving 
a computation of the gravitational wave amplitude, which was shown to be lower than standard 
estimates. Moreover, they showed that for a given mass, the gravitational wave amplitude depends 
only on the deformation of the core. These three-dimensional simulations were made possible 
thanks to the use of spectral methods, particularly for the solution of the Poisson equation for the 
gravitational potential. 

Shock waves give thus difficulties to spectral codes and have been either smoothed with spectral 
vanishing viscosity |112] . or ignored with the code stopping before their appearance. Another idea 
developed first between the Mcudon and Valencia groups was then to use some more appropriate 
techniques for the simulation of shock waves: namely the High- Resolution Shock-Capturing tech- 
niques, also known as Godunov methods (see Living Reviews by Marti and Miiller [145] . and by 
Font [73 ). On the other hand, one wants to keep the fewer degrees of freedom required by spec- 
tral methods for an accurate-enough description of fields, in particular for the solution of elliptic 
equations or for the representation of more regular fields, like the gravitational one. Indeed, even 
in the case where a hydrodynamic shock is present, since it only appears as a source for the metric 
in Einstein's equations, the resulting gravitational field is at least and the spectral series do 
converge, although slower than in the smooth case. Moreover, in a multi-domain approach, if the 
shock is contained within only one such domain, it is then necessary to increase resolution in only 
this particular domain and it is still possible to keep the overall number of coefficients lower than 
the number of points for the HRSC scheme. The combination of both types of methods (HRSC 
and spectral) was first achieved in 2000 by Novak and Ibafiez [159] . They studied a spherically 
symmetric core collapse in tensor-scalar theory of gravity, which is a more general theory than 
General Relativity and allows a priori for monopolar gravitational waves. The system of PDEs 
to be solved resembles the General Relativity one, with the addition of a scalar non-linear wave 
equation for the monopolar dynamical degree of freedom. It was solved by spectral methods, 
whereas the relativistic hydrodynamics equations were solved by Godunov techniques. Two grids 
were used, associated to each numerical technique, and interpolations between both were done at 
every time-step. Although strong shocks were present in this simulation, they were sharply re- 
solved with HRSC techniques and the gravitational field represented through spectral methods did 
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not exhibit any Gibbs-like oscillations. Monopolar gravitational waveforms could thus be given. 
In collaboration with the Garching-hydro group, this numerical technique has been extended in 
2005 to three-dimensions, but in the so-called conformally flat approximation of General Relativ- 
ity (see Sections 15.51 and 15. 6p by Dimmelmeier et al. (69). This approach using spectral methods 
for the gravitational field computation is now sometimes referred as "Marriage des Maillages" 
(French for "grid wedding") and is currently employed by the core-collapse code CoCoNuT of 
Dimmelmeier et al. |681 169) to study general relativistic simulations to a proto-neutron star, with 
a microphysical equation of state as well as an approximate description of deleptonization [70] , 

6.1.2 Collapse to a black hole 

The stellar collapse to a black hole has been a widely studied subject, starting with spherically 
symmetric computations: in the case of dust (matter with no pressure), an analytical solution by 
Oppenheimer and Snyder [163] has been found in 1939. Pioneering numerical works by Nakamura 
and Sato [1541 1155] studied the axisymmetric general relativistic collapse to a black hole; Stark 
and Piran [204] gave the gravitational wave emission from such collapse, in the formalism by 
Bardeen and Piran [19] . Fully general relativistic collapse simulations in axisymmetry have also 
been performed by Shibata [191] , and the first three-dimensional calculations of the gravitational- 
wave emission in the collapse of rotating stars to black holes has been done by Baiotti et al. [17] . 
Recently, Stephens et al. |205] have developed an evolution code for the coupled Einstein Maxwell- 
MHD equations, with the application to the collapse to a black hole of a magnetized, differentially 
rotating neutron stars. 

To our knowledge, all studies of the collapse to a black hole which used spectral methods are 
currently restricted to spherical symmetry. However, in this case and contrary to the core-collapse 
scenario, there is a priori no shock wave appearing in the evolution of the system and spectral 
methods are highly accurate also at modeling the hydrodynamics. Thus, assuming spherical sym- 
metry, the equations giving the gravitational field are very simple, first because of the Birkhoff 's 
theorem, which gives the gravitational field outside the star, and then from the absence of any 
dynamical degree of freedom in the gravitational field. For example, when choosing the radial 
(Schwarzschild) gauge and polar slicing, Einstein equations, expressed within 3-1-1 formalism, turn 
into two first-order constraints which are simply solved by integrating with respect to the radial 
coordinate (see [95]). 

In the work by Gourgoulhon [95] . a Chebyshev tau-method is used. The evolution equations 
for the relativistic fluid variables are integrated with a semi-implicit time scheme and a quasi- 
Lagrangian grid: the boundary of the grid is comoving with the surface of the star, but the grid 
points remains the usual Gauss-Lobatto collocation points fSection r2.3.2p . Due to the singularity- 
avoiding gauge choice, the collapsing star ends in the "frozen-star" state, with the collapse of the 
lapse. This induces strong gradients on the metric potentials, but the code is able to follow the 
collapse down to very small values of the lapse, at less than 10~^. The code is very accurate at 
determining whether a star at equilibrium is unstable, by triggering the physical instability from 
numerical noise at very low level. This property has later been used by Gourgoulhon et al. [102] to 
study the stability of equilibrium configurations of neutron stars near the maximal mass, taking 
into account the effect of weak interaction processes. The addition of some inward velocity field to 
initial equilibrium configurations enabled Gourgoulhon [96] to partially answer the question about 
the minimal mass of black holes: can the effective mass-energy potential barrier associated with 
stable equilibrium states be penetrated by stars with substantial inward radial kinetic energy? 
In |96] , Gourgoulhon found the possibility to form a black hole with a starting neutron star which 
was 10% less massive than the usual maximal mass. 

The spectral numerical code developed by Gourgoulhon [95] has been extended to also simulate 
the propagation of neutrinos, coming from thermal effect and non-equilibrium weak interaction 
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processes. With this tool, Gourgoulhon and Haensel |101j have simulated the neutrino bursts 
coming from the collapse of neutron stars, with different equations of state. Another modification 
of this spectral code has been done by Novak [157] . extending the theory of gravity to tensor- 
scalar theories. This allowed for the simulation of monopolar gravitational waves coming from 
the spherically symmetric collapse of a neutron star to a black hole [157) . From a technical point 
of view, the solution of a non-linear wave equation on curved spacetime has been added to the 
numerical model. It uses an implicit second-order Crank-Nicolson scheme for the linear terms 
and an explicit scheme for the non-linear part. In addition, as for the hydrodynamics, the wave 
equation is solved on a grid, partly comoving with the fluid. The evolution of the scalar field 
showed that the collapsing neutron star "expelled" all of its scalar charge before the appearance 
of the black hole. 

6.1.3 Relativistic stellar pulsations 

Oscillations of relativistic stars are often studied as a time-independent, linear eigenvalue prob- 
lem |130j . Nevertheless, numerical approaches via time evolutions have proved to bring interesting 
results, as obtained by Font et al. [78] for the first quasi-radial mode frequencies of rapidly rotating 
stars in full general relativity. Nonlinear evolution of the gravitational radiation driven instability 
in the r-modes of neutron stars has been studied by many authors (for a presentation of the phys- 
ical problem, see Section 13 of [?'). In particular, the first study of nonlinear r-modes in rapidly 
rotating relativistic stars, via three-dimensional general-relativistic hydrodynamic evolutions has 
been done by Stergioulas and Font [206) . Different approaches doing numerical hydrodynamic sim- 
ulations in Newtonian gravity have been performed by Lindblom et al. |139j . with an additional 
braking term, as by Villain and Bonazzola |224j (see hereafter). 

Because of their very high accuracy, spectral methods are able to track dynamical instabilities 
in the evolution of equilibrium neutron star configurations, as shown in the previous section with 
the works of Gourgoulhon et al. [951 1102] . In these works, when the initial data represented a stable 
neutron star, some oscillations appeared, which corresponded to the first fundamental mode of the 
star. As another illustration of the accuracy, let us mention the work by Novak jl56| . who followed 
the dynamical evolution of unstable neutron stars in tensor-scalar theory of gravity. The instability 
is linked with the possibility for these stars to undergo some "spontaneous scalarization" , meaning 
that they could gain a very high scalar charge, whereas the scalar field would be very weak (or 
even null) outside the star. Thus, for a given number of baryons there would be three equilibria 
for a star: two stable ones with high scalar charges (opposite in sign) and an unstable one with a 
weak scalar charge. Starting from this last one, the evolution code described in |157j was able to 
follow the transition to a stable equilibrium, with several hundreds of damped oscillations for the 
star. This damping is due to the emission of monopolar gravitational waves, which carry away the 
star's kinetic energy. The final state corresponds to the equilibrium configuration, independently 
computed by a simple code solving the generalized Tolman-Oppenheimer-Volkoff system with a 
scalar field, up to 1% error, after more than 50,000 time-steps. These studies could be undertaken 
with spectral methods because in these scenarios the flow remains subsonic and one does not expect 
any shock to be formed. 

It is therefore quite natural to try and simulate stellar pulsations using spectral methods. 
Unfortunately, there have been only a few such studies, which are detailed hereafter. The work 
by Lockitch et al. [142] has studied the inertial modes of slowly rotating stars in full general 
relativity. They wrote down perturbation equations in the form of a system of ordinary differential 
equations, thanks to a decomposition onto vector and tensor spherical harmonics. This system is 
then a nonlinear eigenvalue problem for the dimensionless mode frequency in the rotating frame. 
Equilibrium and perturbation variables are then expanded onto a basis of Chebyshev polynomials, 
taking into account the coordinate singularity at the origin and parity requirements. The authors 
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were therefore able to determine the values of the mode frequency making the whole system 
singular and looked for eigenfunctions, through their spectral decomposition. They found that 
inertial modes were slightly stabilized by relativistic effects. 

A different and maybe more natural approach, namely the time integration of the evolution 
equations, has been undertaken by Villain et al. [2241 1225] with a spectral hydro code, in spherical 
coordinates. The code solves the linearized Euler or Navier-Stokes equations, with the anelastic 
approximation. This approximation, which is widely used in other fields of astrophysics and 
atmospheric physics, consists in neglecting acoustic waves by assuming that time derivatives of 
the pressure and the density perturbations are negligible. It allows for a characteristic time which 
is not set by acoustic propagation time, but is much longer and the time-step can be chosen so 
as to follow the inertial modes themselves. In their 2002 paper |224j . the authors study inertial 
modes (i.e. modes whose restoring force is the Coriolis force, among which the r— modes |5j) in 
slowly rotating polytropes with 7 = 2, in the linear regime. First, this is done in the framework of 
Newtonian gravity, where the anelastic approximation implies that the Eulerian perturbations of 
the gravitational potential do not play any role in the velocity perturbations. Second, they study 
the relativistic case, but with the so-called Cowling approximation, meaning again that the metric 
perturbations are discarded. In both regimes and trying different boundary conditions for the 
velocity field at the surface of the star, they note the appearance of a polar part of the mode and 
the "concentration of the motion" near the surface, showing up in less than 15 periods of the linear 
r— mode. A more recent work |225) deals with the study of gravity modes, in addition to inertial 
modes, in neutron stars. The interesting point of this work is the use of quite a realistic equation of 
state for nuclear matter, which is valid even when the beta equilibrium is broken. The authors were 
thus able to show that the coupling between polar and axial modes is increasing with the rotation 
of the star, and that the coupling of inertial modes with gravity modes in non-barotropic stars can 
produce fast energy exchanges between polar and axial parts of the fluid motion. From a numerical 
point of view, one of the key ingredients is the solution of the vector heat equation, coming from 
the Euler or Navier-Stokes equations. This is done by a poloidal-toroidal [47] decomposition of the 
velocity field onto two scalar potentials, which is very natural within spectral methods. Moreover, 
to ensure the correct analytical behavior at the origin, all scalar quantities are projected at each 
time-step to a modified Legendre function basis. 

More recently, a complete non-linear study of rotating star pulsations has been set by Dim- 
melmeier et al. [71]. They used the general- relativistic code CoCoNuT (see above. Section [6. l.ip 
in axial symmetry, with a HRSC hydrodynamic solver, and spectral methods for the simplified 
Einstein equations (conformally flat three- metric). They noted that the conformal flatness condi- 
tion did not have much effect on the dynamics, when comparing with the Cowling approximation. 
Nevertheless, they found that differential rotation was shifting the modes to lower frequencies and 
they confirmed the existence of the mass-shedding induced damping of pulsations. 

6.2 Vacuum and black hole evolutions 

If one wants to simulate the most interesting astrophysical sources of gravitational radiation, one 
must have a code able to follow, in a stable manner, gravitational waves themselves on a background 
spacetime. It has been observed by all numerical relativity groups that the stability of a numerical 
code, which solves Einstein field equations, does not only depend on the numerical algorithm, but 
also on the particular formulation of the equations. Successes in the simulations of binary systems 
of compact objects in General Relativity (see Section 16. 3p are also due to numerous studies and 
advances in the formulations of Einstein equations. The methods known at present that work for 
the numerical evolution of binaries are the generalized harmonic coordinates |821 1861 1177) and the 
so-called BSSN (for Baumgarte-Shapiro-Shibata-Nakamura [551 ll^5j ). In particular, these two 
formulations of the field equations have the important property that constraint violating modes (see 
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discussion hereafter) stay at a reasonable level during the evolution. Generalized harmonic gauge 
needs constraint damping terms, with a particular method suited for harmonic evolution which was 
proposed by Gundlach et al. that enabled Pretorius to evolve black hole spacetimes [17711176] . 

It is therefore a crucial step to devise such a stable formulation, and more particularly with 
spectral methods, because they add very little numerical dissipation and thus, even the smallest 
instability is not dissipated away and can grow up to unacceptable level. The situation becomes 
even more complicated with the setup of an artificial numerical boundary at a finite distance from 
the source, needing appropriate boundary conditions to control the physical wave content, and 
possibly to limit the growth of unstable modes. All these points have been extensively studied 
since 2000 by the Caltech/Cornell groups and their pseudospectral collocation code SpEC |125| 
[T271 nm [T871 [T881 [T201 [TMl [1371 149]; they have been followed in 2004 by the Meudon group [37] 
and in 2006 by Tichy [2T9] . 

Next, it is necessary to be able to evolve black holes. Successful simulation of binary black holes 
have been performed using the so-called black- hole puncture technique [55l [18] . Unfortunately, the 
dynamical part of Einstein fields are not regular at the puncture points and it seems difficult to regu- 
larize them so as to avoid any Gibbs-like phenomenon using spectral methods. Therefore punctures 
are not generally used for spectral implementations; instead the excision technique is employed, 
removing part of the coordinate space inside the apparent horizon. There is no need for boundary 
condition on this new artificial boundary, provided that one uses free-evolution scheme [1881 , solv- 
ing only hyperbolic equations. In the considered scheme, and also for hydrodynamic equations, 
one does not need to impose any boundary condition, nor do any special treatment on the excision 
boundary, contrary to finite difference techniques, where one must construct special one-sided dif- 
ferencing stencils. On the other hand, with a constrained scheme, elliptic-type equations are to be 
solved [37J and, as for initial data (see Sections 15.31 and 15. 6p boundary conditions must be provided 
e.g. on the apparent horizon, from the dynamical horizon formalism [105) . 

Finally, good outer boundary conditions must be devised, which are at the same time mathe- 
matically well-posed, consistent with the constraints and prevent as much as possible reflections 
of outgoing waves. In that respect, quite complete boundary conditions have been obtained by 
Buchman and Sarbach [55] , 

6.2.1 Formulation and boundary conditions 

Several formulations have been proposed in the literature for the numerical solution of Einstein 
equations, using spectral methods. The standard one is the 3-1-1 (a.k.a. Arnowitt-Deser-Misner 
- ADM) formalism of general relativity [HI 1229) (for a comprehensive introduction, see the lec- 
ture notes by Gourgoulhon [H]), which has been reformulated into the BSSN [^ 1195) for better 
stability. But first, let us mention an alternative characteristic approach based on expanding null 
hypersurfaces foliated by metric 2-spheres developed by Bartnik [20] . This formalism allows for a 
simple analysis of the characteristic structure of the equations and uses the standard "edth" (3) op- 
erator on to express angular derivatives. Therefore, Bartnik and Norton [21] used spin-weighted 
spherical harmonics (see Section I3.2.2p to numerically describe metric fields. 

Coming back to the 3-1-1 formalism, Einstein's equations split into two subsets of equations. 
First, the dynamical equations specifying the way the gravitational field evolves from one time-slice 
to the next; then, the constraint equations which must be fulfilled on each time-slice. Still, it is 
well-known that for the Einstein system, as well as for the Maxwell's equations of electromagnetism, 
if the constraints are verified on the initial time-slice, then the dynamical equations guarantee that 
they shall be verified in the future of that time-slice. Unfortunately, when numerically doing such 
free evolution, i.e. solving only for the dynamical equations, small violations of the constraints due 
to round-off errors appear to grow exponentially (for an illustration with spectral methods, see 
e.g. [188[ 1219] ). The opposite strategy is to discard some of the evolution equations, keeping the 
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equations for the two physical dynamical degrees of freedom of the gravitational field, and to solve 
for the four constraint equations: this is a constrained evolution [57] , 

The advantages of the free evolution schemes are that they usually allow for a writing of the 
Einstein's equations in the form of a strongly- or symmetric-hyperbolic system, for which there are 
many mathematical theorems of existence or well-posedness. In addition, it is possible to analyze 
such systems in terms of characteristics, which can give very clear and easy-to-implement boundary 
conditions [126] . Using finite-differences numerical methods, it is also very CPU-time consuming 
to solve for constraint equations, which are of elliptic type, but this is not the case with spectral 
methods. On the other hand, constrained evolution schemes have by definition the advantage 
of not being subject to constraint-violation modes. Besides, the equations describing stationary 
space-times are usually elliptic and are naturally recovered when taking the steady-state limit of 
such schemes. Finally, elliptic PDEs usually do not exhibit instabilities and are known to be well- 
posed. To be more precise, constrained evolution using spectral methods has been implemented 
by the Meudon group [37], within the framework of BSSN formulation. Free-evolution schemes 
have been used by Tichy [219] (with the BSSN formulation) and by the Caltech/Cornell group, 
which has developed their Kidder-Scheel-Teukolsky (KST) scheme |127] and have later used the 
Generalized-Harmonic (GH) scheme [137] . The KST scheme is in fact a 12-parameters family of 
hyperbolic formulations of Einstein's equations, which can be fine-tuned in order to stabilize the 
evolution of e.g. black hole spacetimes [188K 

Even when doing so, constraint-violating modes grow exponentially and basically three ways 
of controlling their growth have been studied by the Caltech/Cornell group. First, the addition 
of multiples of the constraints to the evolution system in a way to minimize this growth. The 
parameters linked with these additional terms are then adjusted to control the evolution of the 
constraint norm. This generalized version of the dynamical constraint control method used by 
Tiglio et al. [221] . has been presented by Lindblom et al. [138] . and tested on a particular rep- 
resentation of the Maxwell equations. Second, the same authors devised constraint preserving 
boundary conditions from those of Calabrese et al. , where the idea was to get maximally dis- 
sipative boundary conditions on the constraint evolution equations [138] 1126] . This second option 
appeared to be more efficient, but still did not completely eliminate the instabilities. Finally, bulk 
constraint violations cannot be controlled by constraint-preserving boundary conditions alone, so 
Hoist et al. |120l derived techniques to project at each time-step the solution of the dynamical 
equations onto the constraint sub-manifold of solutions. This method necessitates the solution of 
a covariant inhomogeneous Helmholtz equation to determine the optimal projection. Nevertheless, 
the most efficient technique seems to be the use of the GH formulation, which also incorporates 
multiples of the constraints thus exponentially suppressing bulk constraint violation, together with 
constraint-preserving boundary conditions 137 . 

Boundary conditions are not only important for the control of the constraint-violation modes 
in free evolutions. Because they cannot be imposed at spatial infinity (see Section l3.1.2p . they 
must be completely transparent to gravitational waves and prevent any physical wave from en- 
tering the computational domain. A first study of interest for numerical relativity has been done 
by Novak and Bonazzola [158] . where gravitational waves are considered in the wave zone, as 
perturbations of flat spacetime. The specificity of gravitational waves is that they start at the 
quadrupole level (£ — 2) in terms of spherical harmonics expansion. Standard radiative boundary 
conditions (known as Sommerfeld boundary conditions |202] ) being accurate only for the £ = 
component, a generalization of these boundary conditions has been done to include quadrupolar 
terms 158] . They strongly rely on the spectral decomposition of the radiative field in terms of 
spherical harmonics and on spherical coordinates. More specific boundary conditions for the Ein- 
stein system, in order to control the influx of the radiative part of the Weyl tensor, have been 
devised by Kidder et al. [126] for the KST formulation, generalizing earlier work by Stewart [208] 
and Calabrese et al. [SI]. They were then adapted to the GH formulation by Lindblom et al. [137] . 
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Rinne [180] has studied the well-posedness of the initial-boundary value problem of the GH formu- 
lation of Einstein equations. He has considered first-order boundary conditions which essentially 
control the incoming gravitational radiation through the incoming fields of the Weyl tensor. He has 
also tested the stability of the whole system with a collocation pseudo-spectral code simulation a 
Minkowski or Schwarzschild perturbed spacetime. Generalizing previous works, a hierarchy of ab- 
sorbing boundary conditions has been introduced by Buchman and Sarbach [53] , which have then 
been implemented in the Caltech/Cornell SpEC code by Ruiz et al. [182] . together with new sets 
of absorbing and constraint-preserving conditions in the generalized harmonic gauge. Ruiz et al. 
have shown that their second-order conditions can control the incoming gravitational radiation, 
up to a certain point. In addition, they have analyzed the well-posedness of the initial-boundary 
value problems arising from these boundary conditions. Rinne et al. [181] have compared various 
methods for treating outer boundary conditions. They have used the SpEC code to estimate the 
reflections caused by the artificial boundary, as well as the constraint violation it can generate. 

6.2.2 Gauges and wave evolution 

The final ingredient before performing a numerical simulation of the dynamical Einstein system is 
the gauge choice. For example, the analytical study of the linearized gravitational wave in vacuum 
has been done with the harmonic gauge, for which the coordinates {x^} verify the scalar covariant 
wave equation 

H^, = g^.V^V'^a;'' = 0. (185) 

This is the definition of the form i/^, where g^i, is the metric and Vo- the associated covariant 
derivative. Recent works by the Caltech/Cornell group used the GH formulation in which the 
gauge choice is achieved through the specification of Hf^ as an arbitrary function of {x^} and 
Qfii,, which can be set for instance to its initial value [189] . Still, it is with the KST formulation, 
and with lapse and shift set from the analytic values, that Boyle et al. [49] have submitted the 
Caltech/Cornell SpEC code to the so-called "Mexico City tests" [1]. These are a series of basic 
numerical relativity code tests to verify their accuracy and stability, including small amplitude 
linear plane wave, gauge wave and Gowdy spacetime evolutions. These tests have been passed by 
the Caltech-Cornell code, using Fourier basis for all three Cartesian coordinates, and a fourth-order 
Runge-Kutta time-stepping scheme. In the particular case of the linear plane wave, they exhibited 
the proper error behavior, as the square of the wave amplitude, because all non-linear terms are 
neglected in this test. The authors have also shown that the use of filtering of the spherical 
harmonics coefficients was very effective in reducing nonlinear aliasing instabilities. Gauge drivers 
for the GH formulation of Einstein equations have been devised by Lindblom et al. [136] . They 
provide an elegant way of imposing gauge conditions that preserve hyperbolicity for many standard 
gauge conditions. These drivers have been tested with the SpEC code. 

Within the constrained formulation of Einstein's equations, the Meudon group has introduced a 
generalization of the Dirac gauge to any type of spatial coordinates [37' . Considering the conformal 
3-1-1 decomposition of Einstein's equations, the Dirac gauge requires that the conformal 3-metric 
7'^ (such that det7*-' — 1) be divergence-free with respect to the flat 3-metric (defined as the 
asymptotic structure of the 3-metric and with the associated covariant derivative D) 

Dif^ = 0. (186) 

Time coordinate is set by the standard maximal slicing condition. These conditions turn to be 
dynamical gauge conditions: the lapse and the shift are determined through the solution of elliptic 
PDEs at each time-step. With this setting, Bonazzola et al. have studied the propagation of a 
three-dimensional gravitational wave, i.e. the solution of the fully nonlinear Einstein equations in 
vacuum. Their multidomain spectral code based on the Lorene library f5S] was able to follow 
the wave using spherical coordinates, including the (coordinate) singular origin, and to let it out 
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of the numerical grid with transparent boundary conditions [158! . Evolution was performed with 
a second-order semi-implicit Crank-Nicolson time scheme, and the elliptic system of constraint 
equations was solved iteratively. Since only two evolution equations were solved (out of six), the 
other were used as error indicators and proved the awaited second-order time convergence. A 
preliminary analysis of the mathematical structure of the evolution part of this formalism done by 
Cordero et al. (6^ has shown that the choice of Dirac's gauge for the spatial coordinates guarantees 
the strongly hyperbolic character of that system as a system of conservation laws. 

6.2.3 Black hole spacetimes 

As stated at the beginning of Section 16.21 the detailed strategy to perform numerical simulations 
of black hole spacetimes depends on the chosen formulation. With the characteristic approach, 
Bartnik and Norton [21] modeled gravitational waves propagating on a black hole spacetime, in 
spherical coordinates but with a null coordinate z = t — r. They interestingly combined a spectral 
decomposition on spin-weighted spherical harmonics for the angular coordinates and an eighth- 
order scheme using spline convolution to calculate derivatives in the r or z direction. Integration 
in these directions was done with a fourth- or eighth-order Runge-Kutta method. For the spectral 
part, they had to use Orszag's 2/3 rule [56] for anti-aliasing. This code achieved 10~^ as global 
accuracy and was able to evolve the black hole spacetime up to 2; = 55 Af. More recently, Tichy 
has evolved a Schwarzschild black hole in Kerr-Schild coordinates in the BSSN formulation, up 
to t ~ lOOM |219j . He used spherical coordinates in a shell-like domain, excising the interior of 
the black hole. The expansion functions are Chebyshev polynomials for the radial direction, and 
Fourier series for the angular ones. 

Most successful simulations in this domain have been performed by the Caltech/Cornell group, 
who seem to be able to stably evolve forever not only a Schwarzschild, but also a Kerr black hole 
perturbed by a gravitational wave pulse |137j , using their GH formulation with constraint damping 
and constraint-preserving boundary conditions. However, several attempts have been reported by 
this group before that, starting with the spherically symmetric evolution of a Schwarzschild black 
hole by Kidder et al. [128] . Problems had arisen when trying three-dimensional simulations of such 
physical systems with the new parameterized KST formalism [1271 . Using spherical coordinates in 
a shell-like domain, the authors decomposed the fields (or Cartesian components for tensor fields) 
on a Chebyshev radial base and scalar spherical harmonics. The integration in time was done 
using a fourth-order Runge-Kutta scheme and the gauge variables were assumed to keep their 
analytical initial values. The evolution was limited by the growth of constraint-violating modes at 
t ^ lOOOM. With a fine-tuning of the parameters of the KST scheme, Scheel et al. [188] have been 
able to extend the lifetime for the numerical simulations to about 8000M. On the other hand, when 
studying the behavior of a dynamical scalar field on a fixed Kerr background, Scheel et al. [187] 
managed to get nice results on the late time decay of this scalar field. They had to eliminate high- 
frequency numerical instabilities, with a filter on the spherical harmonics basis, following again 
Orszag's 2/3 rule [56] and truncating the last third of coefficients. It is interesting to note that 
no filtering was necessary on the radial (Chebyshev) basis functions. Somehow more complicated 
filtering rule has been applied by Kidder et al. |126j , when trying to limit the growth of constraint- 
violation in three-dimensional numerical evolutions of black hole spacetimes, with appropriate 
boundary conditions. They have set to zero the spherical harmonics terms with £ > imax — 3 in 
the tensor spherical harmonics expansion of the dynamical fields. The stable evolutions reported 
by Lindblom et al. |137] thus might be due to the following three ingredients: 

• GH formalism, exponentially suppressing all small short-wavelength constraint violations, 

• constraint-preserving boundary conditions, 

• filtering of spherical harmonics spectral coefficients. 
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6.3 Binary systems 

As seen in Section 16.21 not many groups using spectral methods are able to solve all the three- 
dimensional Einstein equations in a stable way. When dealing with black holes, the situation is 
even worse. Only very recently, the Caltech/Cornell group succeeded in following 16 orbits, merger 
and ring-down of an equal- mass non-spinning binary black hole system 186 . Moreover, we report 
on three recent partial works on the field using spectral methods, dealing with each type of bi- 
nary system (neutron stars and/or black holes) and leave space for future studies on this rapidly 
evolving field. We note, of course, that successful numerical evolutions of such systems have been 
performed with other numerical methods, which we very briefly summarize here. First successful 
fuUy-relativistic simulation of binary neutron stars has been obtained by Shibata et al. |1971 1193j 
and now, more groups are also able to study such systems: the Louisiana State University (LSU) 
group [1] and the Albert Einstein Institute (AEI, Golm) group [IS]. We also mention here the simu- 
lations with more detailed microphysics by Oechslin and Janka [162] , although with the conformally 
flat approximation, and those of Liu et al. |141j evolving magnetized neutron star binaries. Shibata 
and Uryu [198| I199j have successfully evolved black hole-neutron star binaries, using the puncture 
technique for the modeling of the black hole. As far as black hole binary systems are concerned, 
after many years of hard work and codes evolving the binary system for a restricted time, a first 
stable simulation up to the merger phase has been performed by Pretorius |176| who used general 
harmonic coordinates together with constraint-damping terms and a compactification of spatial 
infinity. He also used the excision technique for a region surrounding the singularity inside the 
horizon. This first success was followed a few moths later by the Texas /Brownsville group [55] and 
the NASA/Goddard group |18j . using very different techniques, namely BSSN with moving punc- 
tures and "l-|-log" slicing together with 'T-driver" shift condition. These techniques have rapidly 
become standards for many other groups, which are now able to stably evolve binary black holes, 
as the AEI/LSU collaboration |174l . group in Jena [95, at Pennsylvania State University [114j and 
Florida Atlantic University [220] . The results have reached high-confidence level and it was possi- 
ble to compare gravitational waveforms obtained numerical evolution to post-Newtonian template 
families [166] . 

6.3.1 Binary neutron stars 

Numerical simulations of the final stage of inspiral and merger of binary neutron stars has been 
performed by Faber et al. [75] . who have used spectral methods in spherical coordinates (based on 
LORENE library [99]) to solve the Einstein equations in the conformally flat approximation (see Sec- 
tions [5] and [HITTT]). The hydrodynamic evolution has been computed using a Lagrangian smoothed 
particle hydrodynamics (SPH) code. As for the initial conditions, described in Section 15. 5i the 
equations for the gravitational field reduce, in the case of the conformally flat approximation, to a 
set of five non-linear coupled elliptic (Poisson-type) PDEs. The considered fields (lapse, shift and 
confornial factor) are "split" into two parts, each component being associated to one of the stars 
in the binary. Although this splitting is not unique, the result obtained is independent from it 
because the equations with the complete fields are solved iteratively, for each time-step. Bound- 
ary conditions are imposed to each solution of the field equations at radial infinity, thanks to a 
multidomain decomposition and a u — 1/r compactification in the last domain. The authors used 
^ 10^ SPH particles for each run, with an estimated accuracy level of 1-2%. Most of the CPU 
time was spent in calculating the values of quantities known from their spectral representation, 
at SPH particle positions. Another difficulty has been the determination of the domain bound- 
ary containing each neutron star, avoiding any Gibbs phenomenon. Because the conformally fiat 
approximation discards gravitational waves, the dissipative effects of gravitational radiation back 
reaction were added by hand. The authors used the slow-motion approximation [227] to induce a 
shrinking of the binary systems, and the gravitational waves were calculated with the lowest-order 
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quadrupole formulas. The code has passed many tests and, in particular, they have evolved several 
quasi-equilibrium binary configurations without adding the radiation reaction force with resulting 
orbits that were nearly circular (change in separation lower than 4%. The code was thus able to 
follow irrotational binary neutron stars, including radiation reaction effects, up to the merger and 
the formation of a differentially rotating remnant, which is stable against immediate gravitational 
collapse for reasonably stiff equations of state. All the results agreed pretty well with previous 
relativistic calculations. 

6.3.2 Black hole-neutron star binaries 

A similar combination of numerical techniques has been used by Faber et al. [74] to compute the 
dynamical evolution of merging black hole-neutron star binaries. In addition to the conformally 
flat approximation and similarly to Taniguchi et al. [210] . the authors considered only the case of 
extremely large mass ratio between the black hole and the neutron star, holding thus the black hole 
position fixed and restricting the spectral computational grid to a neighborhood of the neutron 
star. The black hole surrounding metric was thus supposed to keep the form of a Schwarzschild 
black hole in isotropic coordinates. The neutron star was restricted to low compactness (only a 
few percents) in order to have systems that disrupt well outside the last stable orbit. The system 
was considered to be in corotation and, as for binary neutron stars, the gravitational radiation 
reaction was added by hand. As stated above, the numerical methods used SPH discretization to 
treat dynamical evolution of matter, and the spectral library Lorene [99] to solve the Einstein 
field Poisson-like equations in the conformally flat approximation. But here, the spectral domains 
associated with the neutron star did not extend to radial infinity (no compactified domain) and 
approximate boundary conditions were imposed, using multipole expansion of the fields. The main 
reason is that the black hole central singularity could not be well described on the neutron star 
grid. 

The authors have studied the evolution of neutron star-black hole binaries with different poly- 
tropic indices for the neutron star matter equation of state, the initial data being obtained as 
solutions of the conformal thin-sandwich decomposition of Einstein equations. They found that, 
at least for some systems, the mass transfer from the neutron star to the black hole plays an im- 
portant role in the dynamics of the system. For most of these systems, the onset of tidal disruption 
occurred outside the last stable orbit, contrary to what had been previously proposed in analytical 
models. Moreover, they have not found any evidence for significant shocks within the body of the 
neutron star. This star possibly expanded during the mass loss, eventually loosing mass outward 
and inward, provided that it was not too far within the last stable orbit. Although the major part 
of released matter remained bound to the black hole, a significant fraction could be ejected with 
sufficient velocity to become unbound from the binary system. 

6.3.3 Binary black holes 

Encouraging results concerning binary black holes simulations with spectral methods have been 
first obtained by Scheel et al. [189 . They have used two coordinate frames to describe the motion 
of black holes in the spectral grid. Indeed, when using excision technique (punctures are not 
regular enough to be well represented by spectral methods), excision boundaries are fixed on the 
numerical grid. This can cause severe problems when, due to the movement of the black hole, the 
excision surface can become timelike and the whole evolution problem is ill-posed in the absence of 
boundary conditions. So one solution seems to be the use of comoving coordinates, but the authors 
report that the GH formulation they use appear to be unstable with this setting. They therefore 
consider a first system of inertial coordinates (with respect to spatial infinity) to define the tensor 
components in the triad associated with these coordinates; and a second system of comoving (in 
some sense) coordinates. In the case of their binary black hole tests [189) . they define the comoving 
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coordinates dynamically, with a feedback control system that adjusts the moving coordinate frame 
to control the location of each apparent horizon center. 

The spectral code uses 44 domains of different types (spherical and cylindrical shells, rectan- 
gular blocks) to describe the binary system. Most of the numerical strategy to integrate Einstein 
equations is taken from their tests on the GH formulation by Lindblom et al. [137j and have already 
been detailed in Section 16.2.11 the important technical ingredient detailed by Scheel et al. |189j 
is the particular filtering of tensor fields in terms spherical harmonics. The dual-coordinate-frame 
representation can mix the tensor spherical harmonic components of tensors. So, in their filtering 
of the highest-order tensor spherical harmonic coefficients, the authors had to take into account 
this mixing by transforming spatial tensors to a rotating frame tensor spherical harmonic basis 
before filtering and then transforming back to inertial frame basis. This method allowed them to 
evolve binary black hole spacetimes for more than four orbits, until t > GOOMadm- 

However, a central problem has been the capability of the code to follow the merger phase, 
and even though the code was able to compute the inspiral quite accurately, it used to fail just 
before the holes merged. The problem was that, when both black holes were coming close to 
each other, their horizons became extremely distorted and strong gradients would develop in the 
dynamical fields. This has been explained as a gauge effect, coming from the incapacity of the gauge 
condition to react and change the geometry when the two black holes begin to interact strongly, 
and can be seen as a coordinate singularity developing close to the merger. Nevertheless, a cure 
has been found as explained in Scheel et al. |186| . The original gauge is kept until some given time 
and then smoothly changed to a new one, based on the gauge treatment by Pretorius [177} I176j 
(for the lapse): the gauge source function is evolved through a damped, driven wave equation, 
which drives the lapse toward unity and the shift vector toward zero near the horizons. Thus the 
horizons expand in coordinate space and the dynamical fields are smoothed out near the horizons, 
preventing gauge singularities from developing. With this transition of the gauge condition, the 
evolution of the black holes can be tracked until the formation of a common horizon encompassing 
both black holes. Then, the evolution of this single-distorted dynamical black hole is achieved 
by first interpolating all variables onto a new computational domain containing only one excised 
region, then by choosing a new comoving coordinate system, and finally by modifying again the 
gauge condition to drive the shift vector to a time-independent state. 

This new gauge conditions have allowed Scheel et al. [186] to follow the inspiral during 16 orbits, 
then the merger and the ring-down phase of an equal-mass non-spinning binary black hole system. 
They were able to compute the mass and the spin of the final black hole with very high accuracy 
(10~^ and lO"** relative accuracy for the mass and spin respectively), and to extract the physical 
waveform up accurately to 0.01 radian in phase. This is the first spectral numerical simulation of 
the full evolution of a binary black hole system. 
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7 Conclusions 



We would like to conclude our overview of spectral methods in numerical relativity by pointing a 
few items that we feel are particularly interesting. 

7.1 Strengths and weaknesses 

The main advantage of spectral methods, especially with respect to finite difference ones, is the very 
rapid convergence of the numerical approximation to the real function. This implies that very good 
accuracy can usually be reached with only a moderate number of points. This obviously makes the 
codes both faster and less demanding in memory. Various examples of convergence can be found 
in Section [51 However, this rapid convergence is only achieved for C°° functions. Indeed, when 
the functions are less continuous, spurious oscillations appear and the convergence only follows 
a power-law. In the case of discontinuous functions, this is known as the Gibbs phenomenon 
(see the extreme case of Figure [TT|) . Gibbs- like phenomena are very likely to prevent codes from 
converging or to make time evolutions unstable. So spectral methods must rely heavily on domain 
decomposition of space and the domains must be chosen so that the various discontinuities lie 
at the boundaries. Because of this, spectral methods are usually more difficult to implement 
than standard finite differences (see for instance the intricate mappings of [T). The situation is 
even more complicated when the surfaces of discontinuities are not known in advance or have 
complicated shapes. 

Spectral methods are also very efficient at dealing with problems that are related to coordinate 
singularities. Indeed, if the basis functions fulfill the regularity requirements, then all the functions 
will automatically satisfy them. In particular, it makes the use of spherical coordinates much easier 
than with other methods, as explained in Section [3.21 

Another nice feature is the fact that a function can be represented either by its coefficients or 
its values at the collocation points. Depending on the operation one has to perform, it is easier to 
work with the one representation or the other. When working in the coefficients space, one takes 
full advantage of the non-locality of the spectral representation. A lot of operations that would 
be difficult otherwise can then be easily performed, like computing the ratio of two quantities 
vanishing at the same point (see for instance [ST]). 

7.2 Combination with other methods 

Spectral methods have also demonstrated that they can be a valuable tool when combined with 
other methods. For instance, when shocks are present, spectral methods alone have trouble deal- 
ing with discontinuities at the shock interface. However, this can be efficiently dealt with using 
Godunov methods. Such a combination has already been successfully applied to the simulation of 
the oscillations of compact stars in 71 and of core collapse 1165) . 

Spectral methods have also been used in conjunction with a description of the fluid based 
on SPH (smoothed particle hydrodynamics) in the case of binary neutron stars [75] and for the 
merger of one neutron star and one black hole ^74^. In both cases, the fluid is described by an 
ensemble of particles on which forces are applied. Such technique can account for complicated fluid 
configurations, like the spiral arms that are formed during the merger. Such arms would be tricky 
to follow by means of spectral methods alone. On the other hand, the equations describing gravity 
are nicely solved by spectral solvers. 

7.3 Future developments 

Finally, we would like to point out a few of the directions of work that could lead to interesting 
results. Of course, we are not aware of what the other groups have planed for the future. 
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Appropriate choice of coordinates is evidently important. However, for binary systems, rather 
few results have been using the natural choice of the bispherical coordinates. So far, variations 
of such coordinates have only been used by M. Ansorg and collaborators and only in the context 
of initial data |8l IH [7]. We believe that applying those coordinates, or similar coordinates, to 
evolutionary codes could lead to interesting results, in terms of both speed and accuracy. 

The application of spectral methods to theories more complicated the general relativity is 
also something that can be imagined. One of the possible fields of application is the study of 
branes, where there is an additional dimension to spacetime. The fact that spectral methods are 
accurate with relatively few degrees of freedom, makes them a good candidate to simulate systems 
with extra-dimensions. The addition of gauge fields is also something that could be studied with 
spectral methods, to investigate the possibility of "hairy" black holes, for instance. Of course, 
those are just a few leads on what the future applications of spectral methods to the fields of 
relativity might be. 
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